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Abstract 

In  tliis  dissertation,  we  consider  the  solution  of  linear  systems  of  algebraic  equations 
that  arise  from  finite  element  problems.  We  study  some  additive  Schwarz  type  domain 
decomposition  algorithms  for  general,  not  necessarily  selfadjoint,  linear  elliptic  and 
parabolic  PDEs.  We  use  the  GMRES  method  to  solve  the  resulting  linear  systems  of 
equations.  In  each  step  of  the  iteration,  a  number  of  smaller  linear  systems,  which 
correspond  to  the  restriction  of  the  original  problem  to  subregions,  are  solved  instead 
of  the  large  system.  The  number  of  subproblems  can  be  potentially  large  and  these 
methods  are  therefore  promising  for  parallel  computation.  We  estimate  the  rate  of 
convergence  of  the  algorithms.  Numerical  experiments  are  also  reported. 
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Chapter  1 
Introduction 

1.1      An  overview 

In  this  thesis,  we  consider  the  solution  of  linear  systems  of  algebraic  equa- 
tions that  arise  from  finite  element  problems.  We  study  the  additive  Schwarz 
method  and  the  iterative  substructuring  method  for  general,  not  necessar- 
ily selfadjoint,  linear,  second  order,  elliptic  and  parabolic  partial  differential 
equations. 

Domain  decomposition  techniques  are  ver}-  powerful  iterative  methods 
for  such  problems.  In  each  step  of  an  iteration,  a  number  of  smaller  linear 
systems,  which  correspond  to  the  restriction  of  the  original  problem  to 
subregions,  are  solved  instead  of  the  large  original  system.  The  number 
of  subproblems  can  be  potentially  large  emd  these  methods  therefore  are 
promising  for  parallel  computation.  These  are  divide  and  conquer  methods 
and  the  central  mathematical  question  is  to  obtain  a  bound  on  the  rate  of 
convergence  of  the  iteration.  Borrowing  a  term  from  structural  engineering 
computations,  the  subregions  are  often  called  substructures.  We  thus  have 
two  partitions  of  the  region  into  substructures,  which  sometimes  is  used  to 
define  a  coarse,  global  model,  and  into  elements  of  the  finite  element  model. 


The  iterative  methods  most  commonly  used  are  the  conjugate  gradient 
method  for  the  symmetric,  positive  definite  case  and  the  generalized  conju- 
gate residual  methods  (GMRES)  for  the  general,  nonsymmetric  case.  If  the 
symmetric  part  of  the  operator  is  positive  definite,  with  respect  to  a  suitable 
inner  product,  convergence  can  be  guaranteed.  In  this  thesis,  the  rate  of 
convergence  of  all  algorithms  will  be  estimated.  We  show  that  the  additive 
Schwarz  algorithm  is  optimal  for  both  elliptic  and  parabolic  problems  in 
R?  and  Ft?  in  the  sense  that  the  rate  of  convergence  is  independent  of  both 
the  coarse  mesh  size,  defined  by  the  substructures,  and  the  fine  mesh  size. 
The  iterative  substructuring  algorithm  is  not  optimal  in  the  above  sense, 
however,  in  the  B?  case  the  corresponding  rate  of  convergence  depends  only 
mildly  on  the  mesh  parameters.  A  modified  additive  Schwarz  algorithm  is 
also  introduced  for  parabolic  problems  in  B? .  The  rate  of  convergence  is 
independent  of  the  fine  mesh  size. 

So  far,  most  work  on  domain  decomposition  methods  has  concentrated 
on  selfadjoint  elliptic  equations.  Several  optimal  algorithms  axe  known  for 
this  case,  see  [4],  [5],  [2],  [13],  [22],  etc.  For  previous  work  on  nonselfadjoint 
eUiptic  problems,  see  [8],  [20],  [24]. 

The  numerical  schemes  for  the  nonselfadjoint  problems  fall  into  one  of 
the  following  two  categories.  The  first  class  of  methods  consists  of  higher 
order  accurate  methods  such  as  the  Galer kin's  method,  which  is  the  one 
we  will  concentrate  on,  and  of  the  finite  difference  methods  based  on  us- 
ing centered  difference  approximations  for  the  convection  terms.  These 
methods  can  produce  very  good  solutions  in  the  situation  where  either  the 
original  equation  has  a  smooth  solution  or  the  mesh  size  is  relatively  small 
compared  with  the  Reynold's  number  of  the  equation,  roughly  defined  as 
the  ratio  of  the  magnitude  of  the  convection  coefficient  and  the  minimal 


eigenvalue  of  the  diffusion  coefficient  matrix.  A  discussion  can  be  found 
in  [25].  The  disadvantage  of  these  methods  is  that  they  can  produce  large 
oscillation  when  the  solution  is  discontinuous  and  the  mesh  size  is  relatively 
large.  In  the  other  class  we  find  the  streamline  diffusion  methods,  see  [19], 
which  produce  nonoscillating  solutions  when  the  Galerkin's  method  does 
not  work  well.  In  this  thesis,  we  consider  only  the  Galerkin's  method. 

This  thesis  is  organized  as  follows.  In  the  remainder  of  this  chapter, 
we  first  present  some  basic  Sobolev  spaces  and  the  GMRES  method,  which 
is  the  main  iterative  method  we  shall  use  for  solving  our  linear  system  of 
equations.  Then,  we  develop  an  abstract  theory  for  the  additive  Schwarz 
method,  which  unifies  a  number  of  separate  results  developed  for  some 
different  domain  decomposition  algorithms  in  recent  years.  In  Chapter 
2,  we  apply  our  abstract  theory  to  the  additive  Schwarz  method  for  the 
stationary  convection-diffusion  problems  in  R^  and  F{?  and  an  iterative 
substructuring  method  for  the  same  class  of  problems  in  B? .  We  devote 
Chapter  3  to  the  parabolic  convection-diffusion  problems.  Besides  the  two 
algorithms  just  mentioned,  we  also  introduce  a  modified  additive  Schwarz 
method,  which  is  designed  for  parabolic  problems  in  R^ .  Finally,  in  Chapter 
4,  we  present  some  results  of  our  numerical  experiments. 

1.2      Preliminary  materials 
1.2.1      Some  Sobolev  spaces 

We  will  need  a  number  of  Sobolev  spaces.  Let  fi  be  aji  bounded  open  region 
in  R?  or  R^.  By  L^(f2),  we  denote  the  space  of  square  integrable  functions 
on  Q..  (•,  •)  denotes  the  L^(f2)  inner  product.  For  integer  m  >  0,  H"^(Q,)  is 
the  subspace  of  L'^(Q)  for  which, 


||"||//".(n)  =  (^    E    I  {dldxfuix)  \'  dxf"  <  oo.  (1.1) 

\a\<m 

For  5  =  m  +  <7>0,  0<<7<1,  we  define  H'{Q.)  in  terms  of  the  norm, 

Mh'(U)      =      (||«fH'"(Q)+  (1-2) 

JJ^    E   \{{d/dxru{x)  -  {d/dyru{y}){x  -  y)-^'^"^  I'dadyY^'. 

\a\=m 

We  will  also  need  the  seminorm  |  u  \h'{Q)i  which  is  obtained  by  dropping 
all  the  terms  with  derivatives  of  order  less  than  m  in  (1.1)  and  the  first  term 
in  (1.2). 

We  also  use  HqIQ),  s  >  0,  which  is  a  subspace  of  H^{Q)  defined  as  the 
closure  in  H^{Q.)  of  D{Q),  the  space  of  C°°  functions  with  compact  support 
in  17. 

1.2.2      The  GMRES  method 

In  Eisenstat,  Elman  and  Schultz  [14],  the  conjugate  gradient  method  was 
generalized  to  solve  nonsymmetric  linear  systems  of  equations.  The  so 
called  generalized  minimum  residual  method,  GMRES,  has  been  shown  in 
practice  to  be  very  powerful  for  a  large  class  of  problems.  In  their  paper 
the  GMRES  method  and  the  corresponding  theory  in  L^{Q.)  are  given,  but 
in  fact  this  algorithm  and  its  theory  can  be  extended  easily  to  any  Hilbert 
spaces.  We  shall  describe  the  algorithm  and  state  the  theory  without  proof. 
Let  P  be  an  linear  operator  defined  on  the  finite  dimensional  space  i?" 
with  the  inner  product  [•,•],  which  in  practice  is  chosen  to  take  advantage 
of  some  special  properties  of  P.  In  our  applications,  P  is  not  symmetric 
but  positive  definite  with  respect  to  [■,],  i.e.   there  exist  some  x,  y  6  i2", 


such  that 

but  there  exists  a  constant  c  >  0,  such  that 

[x,Px]  >  c[x,x],    \/x  e  R". 

We  are  interested  in  using  the  GMRES  method  to  solve  the  following 
linear  system  of  equations  on  R^: 

Px  =  6, 

where  b  is  given  in  i?". 

The  iteration  begins  with  an  initial  approximate  solution  Xq  G  R"  and 
an  initial  residual  tq  =  b  —  Pxq.  At  the  m""  iteration,  a  correction  vector 
Zm  is  computed  in  the  Krylov  subspace 

fCm{ro)   =  span{ro,Pro,- ■  •  ,P'"~Vo}, 

which  minimizes  the  following  residual 

min      II  6  -  P{xo  +  =)l 

where  the  norm  is  defined  to  be 

II  •  II  =  yn. 

The  m^^  iterate  is  then  Xm  =  xq  +  z^-  It  can  be  shown  that  if  we 
perform  exact  arithmetic  operations,  then  the  solution  would  be  reached  in 
no  more  than  n  iterations. 

GMRES  Algorithm: 


Choose  xo 

Compute  tq   =   b  —  Pxq 

Set  po    =   To 

For  I  =  0  Step  1  Until  Convergence  Do 

a,    =    [r.,Pp,]/[Pp„Pp.] 

^t+l     =     Xi  -\-  QiPi 

n+1    =   r,  -  a,Ppi 

p,+i    =   r,+i  +  ^=0  ^''^Pj 

where  for;  <  i,  if    =    -[Pr.+i,  Pp,]/[Pp„  Pp^]. 

The  work  and  storage  requirements  per  iteration  becomes  very  expen- 
sive if  a  large  number  of  iterations  is  needed  to  achieve  a  given  accuracy. 
Therefore,  to  reduce  the  number  of  iteration  becomes  a  very  important  is- 
sue for  nonsymmetric  problems.  There  are  some  alternative  ways  to  save 
storage,  such  as  the  Orthomin(A;)  and  the  fc-step  restarted  GMRES  method, 
which  will  not  be  considered  in  the  paper. 

According  to  the  theory  of  Eisenstat,  Elman  and  Schultz,  the  rate  of 
convergence  of  the  GMRES  method  can  be  characterized  by  the  ratio  of 
the  minimal  eigenvalue  of  the  symmetric  part  of  the  operator  and  the  norm 
of  the  operator.  Let  us  define  those  two  quantities  as  follows: 

.   Jx,Px] 
Cp  =   mf  — — 

x#0     [x,x] 


and 

r  II  P^  II 

Cp   =   sup 

r#o    II  -r  II 

Let  P*  be  the  adjoint  operator  of  P  with  respect  to  [■,•].  Then  the 
energy  contributed  by  the  symmetric  part  is 

,P  +   P' 

[ x,x\    =    [Px,x\. 

We  have  the  following  theorem  for  the  rate  of  convergence.  In  the  case 
that  [•,  ■]    =   (■,  •),  the  proof  is  given  by  Eisenstat,  Elman  and  Schultz  [14]. 

Theorem  1.1  If  Cp  >  0,  then  the  GMRES  method  converges  and  at  the 
m"'  iteration,  the  residual  is  bounded  as 

||r„||<    (1    -   7^)         II  roll- 

1.3      The  abstract  theory  for  the  additive  Schwarz 
method 

Let  y  be  a  Hilbert  space  consisting  of  real  functions  defined  on  Cl  C  i?"^, 
where  d  is  the  dimension,  with  an  inner  product  {u,v)v  and  the  correspond- 
ing norm  ||u||v-  Let  a;  be  a  subdomain  of  Q.  We  define  (u,i')v(u,)  to  be  the 
inner  product  for  functions  obtained  by  restricting  on  u. 

Suppose  B{-,-)  is  a  bilinear  form  on  V  x  V  and  Ti-)  a,  linear  functional 
on  V  such  that 

(i)  B{-,  •)  is  continuous,  i.e.  there  exists  a  constant  C  >  0  such  that 

\B{u,v)\<C\\u\\v(^)\\v\\vi^^,     Vu,i;6l/, 
where  u  =  {supp  it}  D  {supp  v}. 


(ii)  5(-,  •)  is  y-elliptic,  i.e.  there  exists  a  constant  c  >  0  such  that 

B{u,u)>c\\u\\l,     Vu  G  F. 
(iii)  ^(•)  is  continuous,  i.e.  there  exists  a  constant  C  >  0  such  that 

\:f{u)\<c\\u\\v,   Vu  e  k 

We  define 

Aiu,v)  =  l/2{B(u,v)  +  B{v,u)), 

which  will  be  called  the  symmetric  part  of  ^(•,  •),  and 

Ar{u,v)  =  l/2{B{u,v)  -  B{v,u)), 

which  will  be  called  the  skewsymmetric  part  of  B(-,  ■). 

We  assume,  following  from  the  assumptions  on  B{u,v),  that  A(u,v)  is 
continuous  and  elliptic  in  the  same  sense  as  in  (i)  and  (ii),  which  implies 
that  the  norm  corresponding  to  A{-,-)  is  equivalent  to  the  F-norm.  In  the 
following,  we  shall  use  {•■,•) a  instead  of  (•,  •)v'. 

Our  abstract  variational  equation  reads  as  follows:  Find  u  £  V,  such 
that 

B{u,v)  =  7{v),    \/ve  V.  (1.3) 

In  order  to  define  the  additive  Schwaxz  method,  we  assume  that  there 
exists  a  decomposition  of  V.  Let  Vi,  i  =  0,-  •  • ,  N,  be  subspaces  of  V,  such 
that 

N 
1=0 

Moreover,  we  assume  that  there  exists  a  constant  Co  >  0  such  that  for  any 
V  E  V,  there  exist  Vi  E  Vi,    i  =  0,-  ■  ■  ,N,  such  that 

N 

v  =  J2''i  (1-4) 

i=0 


and 

N 


Lll^.|lv'(..)^^olMI'v,     Vi;ey,  (1.5) 

t=0 

where  a;,  is  the  support  of  V,  defined  as 

(jj,  =  {x\  e  Q,3u  e  V,  such  that  u{x)  ^  0}. 

Note  that  the  constant  Co  may  depend  on  the  number  of  subregions  A''  and 
cdso  some  other  parameters  of  V',  which  may  be  introduced  in  practical 
appHcations.  By  using  this  function  space  decomposition  {Vi}  we  spht  the 
original  problem  defined  on  V,  which  usually  has  a  large  number  of  degree 
of  freedoms,  into  some  smaller  problems  defined  on  Vi,  which  are  chosen  to 
have  relatively  small  number  of  degree  of  freedoms.  These  small  problems 
are  usually  independent  of  each  other  therefore  can  be  solved  in  parallel. 

We  also  assume  that  the  supports  of  {u;,}  form  a  finite  covering  of  Q,  in 
the  sense  that  there  exists  a  constant  C^  >  0,  such  that 

i:ibllv(..)<C.||t'||^,    yvev.  (1.6) 

1=0 

In  fact  this  constant  C^j  is  the  maximum  number  of  u^s  to  which  a  point 
X  G  f2  can  belong. 

For  each  subspace  V^,  0  <  i  <  A^,  we  define  a  projection 

pB    =    pB.     y    ,  y 

with  respect  to  the  bilinear  form  B{-,  •):  For  any  given  u  ^  V,  P^u  G  V^  is 
the  solution  of  the  problem 

B{p!^u,v)  =  B{u,v),  Vvev,. 

We  introduce  a  mapping  P^  :  V  — >  V  as 

9 


which  is  the  main  operator  that  we  shall  study  in  this  thesis. 

Let  u  be  the  solution  of  (1.3),  and  denote  the  image  of  u  under  the 
mapping  P'^  as 

N 

b  =  ZpS- 

1=0 

It  is  easy  to  see  that  b  can  be  computed  without  knowing  the  solution  u 
itself.  We  compute  Pj^u  by  solving  the  linear  system  of  equations 

for  each  i.  b  is  thus  obtained  by  taking  the  sum  of  the  solutions. 

We  now  introduce  a  new  linear  system  of  equations:  Find  u  6  F,  such 
that 

P^u   =   b.  (1.7) 

We  shall  call  this  equation  the  derived  equation  with  respect  to  the  bilinear 
form  B{-,-)  and  the  decomposition  [Vt]. 

The  following  theorem  can  be  established  trivially. 

Theorem  1.2  If  P'^  is  invertiahle,  then  the  equations  (1.3)  and  (1.7)  have 
the  same  solution. 

The  additive  Schwarz  algorithm  can  be  stated  as 

Additive  Schwarz  Algorithm:  Find  the  solution  u  of  equation  (1.3) 
by  solving  equation  (1.7). 

In  practice,  the  operator  P  usually  corresponds  to  the  sum  of  products 
of  a  large,  sparse  matrix  and  the  inverses  of  some  other  sparse  matrices. 
The  explicit  matrix  for  P^  is  not  known  except  in  some  special  cases,  how- 
ever, the  matrix  vector  multiply  P^u  can  be  computed  by  solving  one 
linear  system  of  equations  for  each  subregion.  Therefore  iterative  methods 


10 


are  natural  candidates  for  problem  (1.7).    The  rate  of  convergence  of  any 
iterative  methods  depends  on  the  conditioning  of  the  operator  P^. 

Theorem  1.3   (1)  There  exists  a  constant  C  >  0,  such  that 

\\p'^u\\a  <  cc^uWa,    VueF. 

(2)  There  exists  a  constant  C\  >  0,  independent  of  Co,  such  that 

\\P''u\u>c,Co'\\u\u    yuev. 

(S)  If  there  exists  0  <  6  <  I,  such  that  \  jV  {u ,  P'^  u)  \<  6B{u,P'^u),    Vu  G 
V,  then 

(u^P^uU  >  C2(l  -  S)Co\u,uU,    Vu  G  V, 

where  C2  >  0  is  independent  of  Co  and  6. 

Proof: 

(1)  We  first  show  that  there  exists  a  constant  C  >  0,  such  that 

N 


E  II  Pf^  IU(..)    <   CCJ  u  \\\,     Vu  G  V.  (1.8) 


By  the  definition  of  the  A-novm, 


N  N 

\A{uj,) 


Ell^^"lli(..)  =  E^(^.V^^),  "^^^v. 

t=0  :=0 


By  the  definition  of  the  projection  P^,  we  have 

B{P^u,PS)  =  B{u,Pfu). 

Since  P^u  is  zero  outside  cj,,  we  have,  by  the  boundedness  assumption  (i), 
that 

11 


Combining  the  above  results  and  applying  Cauchy's  inequality,  we  obtain 


EWPFuwi,,^,  <  c  j:\\-\\Ai.,)-\i: 

t=0  N  <=0  \  1=0 


If  we  cancel  the  common  factor  on  both  sides,  and  use  the  assumption  (1.6), 
we  obtain  the  estimate  (1.8). 

It  now  suffices  to  bound  ||  P'-^u  ||^  by  the  left  hand  side  of  estimate  (1.8). 
Using  the  definition  of  P^,  we  have 

llP^ull^   =  5(PVP^)  =  ZBiP^u^P^'u). 


:=0 


Note  that  the  function  Pj^u  has  its  support  only  in  u;,.  Thus 
BiP'u.Pl'u)  <  C||  P^u  11^,^,)  .  II  PM|^(^_). 
Applying  Cauchy's  inequality  and  the  assumption  (1.6),  we  obtmn 


N 


p'u  wi,  <  cc.i: II  p.^u"' 


1=0 


A{u,,)- 


Combined  with  (1.8),  we  can  conclude  that 


P^u  11^    <   CCJI  u  11^, 


which  proves  the  upper  bound  part  of  the  theorem. 

(2)  We  now  turn  to  the  proof  of  the  lower  bound  part  of  the  Theorem  1.3. 
By  the  V-function  decomposition  cissumption  (1.4),  any  u  6  V  can  be 
written  as 
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N 
u     =    ^U„ 
1=0 

where  u,  G  Vi.   This  decomposition  is  chosen  so  that  it  satisfies  (1.5).   We 
have 

N 


1=0 

Using  the  definition  of  the  projection  P,^,  we  obtain 


1=1 


<  c'i:iiPMu^.)-iit..iu^.,. 

1=0 

Applying  Cauchy's  inequcJity  and  the  bounded  decomposition  assumption 
(1.5),  we  obtain 

liHii  <  c'-Co^f;iip.^u|ii(^.),  (1.9) 

1=0 

where  C  >  0  is  independent  of  Cq.  By  the  definition  of  the  yl-norm, 

»=0  1=0 


=  j:5(u,i^^u). 

1=0 

Because  of  the  Uneaxity  of  5(-,  •),  we  have 

Ell  PF^  nil.,)    =   B{u,P^u)    <   C\\  u  U\  P'u  11^.  (1.10) 

1=0 
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Combined  with  the  estimate  (1.9),  we  have 

ilHL  <  cclwr^uW^, 

where  C  is  independent  of  Co- 
(3)  It  is  easy  to  verify  that 

(u,P^u)^    =   B{u,P'^u)-M{u,P'^u). 

By  the  assumption  on  the  biHnear  form  M{-.  •),  the  right  hand  side  can  be 
bounded  from  below  by 

(l-(5)5(u,P^u). 

Since  P^  is  the  sum  of  the  P^  s,  the  above  expression  equals 

(i-6)f:5(p^,pM 

1=0 

By  assumption  (ii),  it  can  be  bounded  from  below  by 

N 


c{i  - S)Y.\\P>\W.,r 


,iP^-'^ 

t=0 

The  constant  C2(Co)  =  ^(l  —  S)Cq^   >  0  can  be  obtained  by  using  the 
estimate  (1.9)  and  C  >  0  is  independent  of  Cq.  □ 
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Chapter  2 

Methods  for  stationary 
convection-diffusion  problems 


In  this  chapter,  we  study  two  Schwarz  type  methods  for  the  stationary 
convection-diffusion  equations;  the  additive  Schwarz  method  (A^I)  and  the 
iterative  substructuring  method  {lo/^I).  Both  methods  fall  in  the  framework 
provided  by  the  abstract  additive  Schwarz  theory,  however,  the  analysis 
and  computational  implementation  are  different.  We  prove  that  A^idis  an 
optimal  algorithm  in  the  sense  that  the  rate  of  convergence  is  independent 
of  the  mesh  parameters  provided  that  the  size  of  the  substructures  satisfies 
some  conditions.  This  analysis  is  valid  in  both  R^  and  R^.  We  also  give  a 
convergence  rate  analysis  for  the  7s.'\ffor  the  R'^  case.  The  bound  depends 
only  mildly  on  the  mesh  parameters. 

2.1      Introduction 

2.1.1      A  stationary  convection-diffusion  problem 

In  this  section,  we  present  a  continuous  second  order  linear  elliptic  par- 
tial differential  equation,  which  is  not  necessarily  selfadjoint,  in  a  bounded 
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region  in  R^  or  R^.  We  also  give  the  variational  formulation,  basic  as- 
sumptions and  some  classical  regularity  results,  which  are  important  for 
the  development  of  our  theory. 

Let  Q  be  an  open  bounded  polygon  in  R^  or  i?"',  with  boundary  dCl. 
d  denotes  the  dimension  of  the  space,  d  ■=  2  or  d  =  Z.  Consider  the 
homogenous  Dirichlet  boundary  value  problem 

f  Lu     =     /    in     Q 

\      u     =     0     on     aQ,  ^       ' 

where  L  is  a  strongly  elliptic  operator  of  the  following  form 

Lu{x)  =  -  2^  ^~(°«i(2;)-^ )  +  2^  o,{x)— +  c{x)u{x), 

where  a^jix)  =  a_,,(x)  for  all  i,  j  and  x  E  Q- 

Since  the  adjoint  problem  of  (2.1)  will  be  used  in  later  chapters,  we 
shall  write  down  an  integral  formula.  The  adjoint  operator,  denoted  as  L', 
satisfies  the  relation: 

{Lu,v)   =  {u,L'v),     yu,veH^in). 

By  applying  Green's  identity,  we  can  easily  find  the  formula  for  L*, 

L  v  =  -  }_^  ■^-{aijix}-^)  -  2^  -^{b,{x)v)  +  c{x)v. 
.,j^i  9xi  dxj        fr{  ox. 

We  assume  that  /  6  L^{Q.).  The  existence  and  uniqueness  of  the  solu- 
tion of  equation  (2.1),  as  well  as  its  adjoint  equation,  are  well  understood. 
It  is  well  known  that  if  the  domain  boundary  is  not  smooth  enough,  for 
example  on  a  polygonal  domain  and  the  domain  is  not  convex,  we  cannot 
expect,  in  general,  the  solution  to  be  in  H^(Q).  According  to  the  classical 
elliptic  regularity  theory  on  Lipschitz  region,  see  [1]  we  make  our  regularity 

16 


assumption  as  follows:  There  exists  a  constant  7  €  (0, 1/2),  which  depends 
on  the  geometry  of  Q,,  such  that  the  adjoint  equation  has  a  unique  solution 
u  6  H^'^'^{Cl)r\HQ{Ct),  and  furthermore  there  exists  a  constant  C,  such  that 

Mn^-Hn)  <  C||L*u|U.(n).  (2.2) 

The  weak  form  of  equation  (2.1)  reads  as  follows;  Find  u  such  that 

B{u,v)  =  F(v),       VrG^^(Q),  (2.3) 

where  the  bilinear  form  is  defined  by 

and  the  linear  functional  F  is  defined  by 

F{v)   =     f   fvdx   with   V  e  H^iO.). 

We  assume  that  this  bilinear  form  is  bounded  and  strongly  elliptic,  i.e. 
there  exist  two  constants  C  >  0  and  c  >  0,  such  that 

I  5(u,i;)l<Ci|u  11^,(^)11  Hl//:(n).     V^i.t,  e  ^^(fi),  (2.4) 

and 

5(u,u)>c||u  ||^,(f^),     \JueHliQ.).  (2.5) 

If  the  coefficients  of  the  first  order  terms  in  equation  (2.1)  are  identi- 
cally zero,  then  the  differential  operator  L  is  selfadjoint.  To  separate  the 
symmetric  part  and  the  skewsymmetric  part  of  the  bilinear  form  B{-,-),  we 
introduce  two  bilinear  forms 

(1)  The  symmetric  part 
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A{u,v)    =     ^{Biu,v}  +  B(v,u)) 

=    Yltj=i  iQaijdu/dXfdv/dxjdx  +  J^c  —  ^div{b)uvdx 
(2)  The  skewsymmetric  part 

N{u,v)    =     l{B{u,v)-B{v,u)) 

—     IIi'=i /n  ^:<5u/5x,?;cfx   +   ^  fQdiv{b)uvdx. 

From  the  boundness  assumption  on  B{-,-)  and  the  Schwarz  inequaHty, 
it  follows  that  there  exists  a  constant  C  >  0  such  that 

I  A{u,v)  \<  C\\  u  11^,(0)11  ^'  ll//,'(0)'    ^"'^  ^  ^o(i^), 

and  from  the  ellipticity  assumption,  there  exists  another  constant  c  >  0 
such  that 

Remark:  The  boundness  and  ellipticity  of  -4(-,-)  imply  that  the  A- 
norm,  defined  by  JA(-,-),  is  equivalent  to  the  Hq(Q,)  norm. 

We  assume  that  the  bilinear  form  A'^(-,  •)  is  bounded,  i.e.  there  exists  a 
constant  C,  such  that 

I  Niu.v)  \<C\\u  y.(a)\\  V  IIl^(Q),       Vu,t;  G  H',{n). 

This  bound  can  be  established  if  the  coefficients  6,(x)  and  div{b)  are  bounded. 
We  note  that  the  estimates  for  A(-,  •)  and  A^(-,  •)  axe  different.   If  we  look 
at  the  integral  formulas  for  A{-,-)  and  iV(-,-),  it  is  clear  that  the  terms  in 
iV(-,-)  are  one  order  lower  than  the  terms  in  A{-,-).    This  is  a  key  factor 
that  makes  our  proofs  work. 
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2.1.2      The  Galerkin  finite  element  method 

We  solve  equation  (2.3)  by  a  conformal  Galerkin  finite  element  method. 
For  simplicity,  we  use  piecewise  linear  triangular  element  in  R^  and  the 
corresponding  tetrahedronal  element  in  R'^.  In  this  subsection,  we  first 
introduce  a  two  level  triangulation  of  ft  and  the  corresponding  finite  element 
spaces.  We  also  discuss  some  basic  properties  of  the  finite  element  spaces. 

1.  A  Two  level  triangulation 

For  a  given  polygonal  region  Q  £  R'^,  in  the  first  step,  we  define  {fi,} 
to  be  a  regular  finite  element  triangulation  of  Q  where  {f2,}  is  a  set  of  non- 
overlapping  (i-dimensional  simplices,  i.e.  triangles  if  <i  =  2  and  tetrahedra 
if  (f  =  3,  such  that  no  vertex  of  one  simplex  lies  on  an  edge  or  a  face  of 
another  and 

_         N  

fi  =  U  Q,. 
1=1 

Here  N  is  the  number  of  simplices. 

Let  Hi  denote  the  diameter  of.  17,  and  H,  denote  the  diameter  of  the 
largest  inscribed  ball  in  17,.  We  assume  that  the  ratio  Hi/H,  is  uniformly 
bounded  from  above,  i.e.  the  discretization  is  shape  regular. 

We  introduce  the  mesh  parameter 

H   =  maxfifi,  •  •  •  ,.H';v}. 

We  caJl  17,  a  substructure  and  {17,}  the  coarse  mesh  or  .H^-level  subdivision 
of  17. 

In  our  second  step,  we  further  divide  each  substructure  17^  into  smaller 
simplices,  denoted  as  r/,  j  =  I,---  .  We  assume  that  {r/}  form  a  shape 
regular  finite  element  triangulation  in  the  same  sense  as  above.  Let  h''-  be 


19 


the  diameter  of  r/ ,  we  introduce  the  fine  mesh  parameter 

h   =   max(/?^). 

We  call  U,,_,r/  the  fine  mesh  or  /i-level  subdivision  of  fl  . 
Next,  we  define  the  piecewise  hnear  finite  element  function  spaces  over 
both  i7-level  and  /i-level  subdivision  of  Q. 

V     =  {v     I    continuous  on  fi,  u     |n,    hnear  on  fi^,  v"  =  0  on  dO,} 

V    =  {v    I    continuous  on  Q,  y    \^j    hnear  on  r/,  u''  =  0  on  dQ} 
It  is  obvious  that  V^  C  V^.  To  simphfy  the  notations,  we  denote 
A     =    {x  |G  interior  nodes  oi  h  —  level  subdivision} 
A      =    {x  |6  interior  nodes  of  if  —  level  subdivision}. 

Associated  with  each  node  x,  G  A'',  there  is  a  nodal  base  function, 
denoted  as  <?!.f,  such  that  (f>'l{x)  E  V"  and  <f>iixj)  =  6,j.  The  set  {(f>'^}  forms 
the  nodal  base  o{  V^. 

2.  Galerkin  finite  element  approximation 

The  Galerkin  approximation  of  equation  (2.3)  reads  as  follows:  Find 
u''  eV",  such  that 

B{u\v'')  =  F{v''),    yv'^eV^.  (2.6) 

The  existence  and  uniqueness  of  u''  has  been  extensively  studied  in  the 
literature,  see  [1].  By  using  the  nodal  base  functions,  equation  (2.6)  can 
be  transformed  into  a  linear  system  of  equations,  which  is  usually  large 
and  sparse.  It  is  well  known  that  the  efficiency  of  any  iterative  methods 
to  be  used  to  solve  this  linear  system  of  equations  depends  strongly  on  the 
conditioning  of  the  stiffness  matrix  A'/i,  where  Kh  =  {B{(f>^,  ^j))- 
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2.2      Additive  Schwarz  method  for  the  prob- 
lems in  R~  and  R^ 

2.2.1      An  algorithm  and  the  main  results 

The  additive  variant  of  the  Schwarz  alternating  method  was  originally  pro- 
posed by  M.  Dryja  and  O.  Widlund  [11]  and  S.  Nepomnyaschikh  [26]  for 
selfadjoint  stationary  elliptic  problems.  In  this  section,  we  generalize  this 
method  to  the  nonselfadjoint  elliptic  case. 

We  first  form  a  basic  decomposition  of  the  domain  ft  and  we  then  define 
the  projections  which  will  lead  to  our  additive  Schwarz  algorithm. 
1.  Basic  domain  decomposition  and  projections 
In  the  previous  section,  we  introduced  the  H-leve\  subdivision  {fli}  of 
fl.  Since  the  Schwarz  type  domain  decomposition  methods  use  overlapping 
subregions,  we  extend  each  subregion  Q,  to  a  larger  region  Q.[  such  that 

Moreover,  we  assume  that  there  exists  a  constant  a  >  0  such  that 
distaiice(5r2;  fl  Q,dQ,  D  Omega)  >   aH,,    Vi. 

To  simplify  the  notation,  we  denote  Q.q   =   Q. 

We  suppose  that  dQ'^  does  not  cut  through  any  h-level  elements.  We 
make  the  same  constructions  for  the  subregions  that  meet  the  boundary 
except  that  we  cut  off  the  parts  that  are  outside  Q,. 

The  assumption  on  the  existence  of  such  a  constant  a  is  important. 
From  the  estimates,  which  will  be  given  later  in  this  chapter,  we  can  see 
that  the  larger  the  a  the  better  our  lower  bound  on  the  spectrum  will  be. 
But  if  we  increase  the  overlap,  the  size  of  the  subproblems  also  increases, 
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therefore,  the  cost  for  solving  the  subproblems  in  each  iteration  will  be 
increased.  To  balance  the  total  number  of  iteration  and  the  cost  of  solving 
the  subproblems  is  an  important  practical  issue. 

For  each  Q-,  there  is  a  regular  finite  element  subdivision  which  is  nat- 
urally induced  from  the  /i-level  subdivision  of  Q.  The  corresponding  finite 
element  function  space,  denoted  by  V-'',  is  defined  as 

which  can  be  regarded  as  a  subspace  of  V''  if  we  extend  each  function  by 
zero  to  the  complement  of  Q^.  It  is  known  that  this  extension  is  continuous. 
In  the  case  when  5Q,  intersects  5Q,  we  use  the  original  boundary  condition 
on  dQ.[f]Ct. 

To  make  our  notations  simple,  we  let 

When  convenient,  we  can  regard  any  element  of  Vq  as  a  coarse  mesh  inter- 
polation Ihu^  of  some  elements  u'^  G  V^,  which  uses  the  function  values  at 
the  nodes  of  the  coarse  mesh  only. 

It  can  be  easily  seen  that  our  finite  element  function  space  V'^  can  be 
represented  as  the  sum  of  the  TV  +  1  subspaces,  i.e. 

v'^  =  Vo^  +  Vl'  +■■■+  F^. 

Let  Pyh   =   P-^  denote  the  projection  from  the  finite  element  space  V'^ 

i 

to  the  subspace  V^''  with  respect  to  the  bilinear  form  5(-,  •),  and  define  the 
mapping  P^  :   V^  -*  V^  as 

pB    =    pB    +    pB    +...    +    pB 
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This  is  the  main  operator  in  our  further  studies  in  the  chapter. 

It  is  easy  to  see  that  the  computation  of  the  projection  of  an  arbitrary 
function  v^  E  V'  into  the  subspace  V^''  involves  only  the  solution  of  a 
standard  finite  element  linear  system  of  algebraic  equations  on  Q,,  which 
is  a  small  subregion  in  the  case  i  ^  0.  If  i  =  0,  it  is  the  standard  finite 
element  equation  on  the  iJ-level  coarse  space. 

Let  us  denote 

b  =  p^u''  =  Y.pS''- 

1=0 

By  taking  b  as  the  right  hand  side,  we  can  define  our  derived  linear  system 
with  respect  to  the  bilinear  form  B{-,-)  and  the  decomposition  {Vj}  as 
follows, 

P^u''   =  h.  (2.7) 

By  Theorem  1.2,  the  equation  (2.7)  is  equivalent  to  the  Galerkin  equa- 
tion (2.6).  If  the  operator  L  is  not  selfadjoint,  then  P^  is  not  symmetric. 

2.  An  algorithm  and  the  main  theorem 

Due  to  the  special  properties  of  P^,  any  numerical  algorithms  designed 
to  solve  such  a  linear  system  of  equations  must  be  robust  with  respect  to 
asymmetry  and  require  matrix  vector  multiplications  only.  The  GMRES 
method  described  in  Chapter  1  is  one  such  algorithm  which  we  shall  use 
in  this  paper.  Chebyshev  iterative  method  (see  [23])  would  provide  an 
alternative,  which  will  not  be  considered  here. 

Additive  Schwarz  algorithm:  Solving  equation  (2.7)  by  GM- 
RES method  with  the  inner  product  [•,  •]   =   (•,)^. 

According  to  Theorem  1.1,  the  rate  of  convergence  of  this  algorithm  can 

23 


be  estimated  by  certain  spectral  bounds  for  the  operator  P^ .  One  of  the 
main  theorems  of  this  chapter  shows  that  the  operator  P^  of  the  -equation 
(2.7)  is  uniformly  well  conditioned  in  the  sense  that  all  the  bounds  are 
independent  of  the  mesh  parameters  H  and  h.  The  proof  is  given  in  the 
next  section. 

We  adopt  the  following  convention.  All  the  constants,  denoted  as  C,  c,  Cp.  Cp 
etc.,  are  positive  and  independent  of  the  mesh  parameters  H  and  h. 

Theorem  2.1  The  operator  P^  is  uniformly  well  conditioned  in  the  fol- 
lowing sense: 

(1)  There  exists  a  constant  Cp,  such  that 

II  P^u'  lU  <  Cp  II  u'  lU,  Vu''  e  v\ 

(2)  There  exists  a  constant  c,  such  that 

\\P^u'\\a>c\\u'\\a,    yu'eV\ 

(S)  There  exists  a  constant  Hq  >  0,  independent  of  H  and  h,  and  a 
constant  Cp{Ho),  such  that  for  H  <  Hq, 

iu\P^u'')A    >   Cpiu\u^)A,    Vu''gF\ 

Remarks:  (a)  The  operator  Pq  is  very  important.  It  provides  for 
global  information  transportation.  All  the  other  Pfs  are  local  mappings. 
Without  using  P(f ,  in  each  iteration  the  informaiion  travels  only  from  one 
substructure  to  its  next  neighbors.  Therefore,  it  takes  at  least  0(1 /H) 
iterations  for  the  information  to  move  across  the  region. 

(b)  In  the  case  that  B{-,-)  is  symmetric,  this  additive  Schwaxz  algorithm 
is  identical  to  the  one  proposed  by  M.  Dryja  and  O.  Widlund.  Theorem  2.1, 
part  (1)  and  (2),  gives  the  convergence  rate  estimates. 
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(c)  Theorem  2.1  part  (3)  shows  that  if  the  coarse  mesh  size  is  fine 
enough,  then  the  symmetric  part  is  positive  definite,  which  guarantees  that 
the  GMRES  method  converges.  Since  both  constants  Cp  and  Cp  are  inde- 
pendent of  the  mesh  parameters  H  and  h,  then  according  to  Theorem  1.1, 
the  convergence  rate  of  the  GMRES  method  will  not  dependent  on  the  size 
of  the  discrete  problems. 

(d)  The  constant  Hq  determines  the  size  of  the  coarse  mesh  problem. 
Ho  depends  on  the  problem.  In  general,  Hq  decreases  if  we  increase  the 
coefficients  of  the  nonselfadjoint  terms,  while  it  increases  if  we  use  larger 
overlap.  It  also  depends  on  the  shape  of  the  domain  fl.  We  do  not  have 
an  explicit  formula  for  Hq.  However,  we  will  give  some  idea  about  how 
to  determine  this  constant  numerically  later  in  the  chapter  of  numerical 
results. 

(e)  If  the  skewsymmetric  part  of  the  elliptic  operator  vanishes,  then 

Ho   =   +00, 
and  hence  we  have  no  restrictions  on  the  coarse  mesh  size  H. 

2.2.2      The  condition  number  estimates 

In  order  to  prove  Theorem  2.1,  we  need  only  to  show  that  all  the  assump- 
tions for  Theorem  1.3  hold.  Since  {fi,}  is  a  finite  covering  of  Q.,  the  maxi- 
mvmi  number  of  fi.s  to  which  any  point  in  f2  can  belong  is  a  finite  number, 
denoted  as  C^.  It  is  easy  to  verify  that 

f:\\^'\\Mn:)<CJu'\\l,    yu'eV\  (2.8) 

«=o 

i.e.  assumption  (1.6)  holds. 
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Note  that  if  we  increase  the  size  of  the  overlap  between  the  basic  sub- 
region  r2,s,  then  this  constant  C^  increases. 

Next,  we  show  that  the  assumption  (1.5)  holds.  Let  us  introduce  the 
coaxse  mesh  L^  projection  as  follows:  For  any  u^  G  V^,  the  L^  projection 
Ihu^  6  V^  is  the  solution  of  the  following  linear  system 

Some  useful  properties  about  the  L^  projection  are  summarized  in  the 
following  lemma,  which  holds  in  both  R^  and  R'^. 

Lemma  2.1    There  exist  two  constants  Ci  >  0  and  Cj  >  0,  independent  of 
H  and  h,  such  that 

||///u'IU<Ci||u''|U,  Vu''eT^\ 

Proof:  See  [27].  D 

Lemma  2.2   Fot  any  u'^  6  V^,  there  exist  u'^  G  V/^{Q,[),     such  that 

N 
t=0 

Moreover,  there  exists  a  constant  Co,  such  that 

E  II  «f  ll'^(n;)   <  Co^ll  "'  C 

t=0 

Proof:    Following  M.  Dryja  and  O.  Widlund,  we  first  construct  the 
decomposition  and  then  prove  that  it  satisfies  the  estimate. 
For  a  given  u''  G  V*,  we  take  Uq  to  be  the  L^  projection 
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Let  us  denote 

Following  the  paper  of  M.  Dryja  and  O.  Widlund  [11],  the  functions  u'^ 
mentioned  in  Lemma  2.2  can  be  constructed  easily.  We  take 

where  {9,}  defines  a  partition  of  unity  of  Q  and  belongs  to  C^{Q.[)  and  Ik 
is  the  interpolation  operator  at  the  h-level  nodal  points.  We  can  arrange  so 
that  V^t  is  bounded  by  const/ Hi.  By  using  the  linearity  of  J/j,  we  can  easily 
show  that  we  have  obtained  a  correct  decomposition  o{  u^.  According  to 
the  theory  of  M.  Dryja  and  O.  Widlund  [11],  we  have 

which  holds  for  both  2D  and  3D.  We  can  add  up  them  for  i  =  1,  •  •  ■  ,iV 
and  obtain  an  estimate  on  f2. 

According  to  Lemma  2.1,  we  have 

II  ^'  WIhu)    <   CH'\  u"  1^. 
The  proof  follows  by  combining  these  estimates.  □ 

In  order  to  prove  part  (3),   we  need  the  following  estimates  for  the 
skewsymmetric  part. 

Lemma  2.3   There  exists  a  constant  C ,  such  that 

(a)     I  N{P^u\u^  -  P^u^)  I  <   CH-^  II  P^u^  lUII  u"  |U,    Vu"  G  V\ 
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(6)     I  N{P^u\u'  -  Pfu')  \<   CH\\  PS'  IL(«;)||  u'  ||^(^;),    yu'  6  V 

fori  =  !,■■■  ,N. 

Proof:   (a)  By  the  boundness  assumption  on  N(-,-),  we  need  only  to 
show  that 

II  u'  -  P,^u'  U^u)  <   CH-^  II  u'  lU, 

where  C  >  0  does  not  dependent  on  H  and  h.  From  the  identity 

we  can  easily  obtain  the  -4-norm  estimate,  i.e.  there  exists  a  constant  C  >  0 
such  that 

II  u''   -   Po^u'^  lU  <    C  II  k''  IU  . 

The  L^  estimate  is  obtained  by  using  the  so  called  Nitsche  trick.  We  have 

II  Po'u'   -   u"  |U.(.)  =      sup     (Po-'-u\v)^  ^2^^ 

IIHlL2?to       II  V  Wmn) 

For  any  fixed  v  G  L'^{fl),  we  introduce  the  following  auxiliary  variational 
problem:  Find  w,  such  that 

Bicf>,w)   =  (t;,<^),      Wcf>eH^in). 

By  our  assumptions  the  solution  exists  and  satisfies 

II  w  \\m+-,(n)  <  C"  II  i^  Wl^o)  ■ 

If  we  take  the  test  function  <f)  in  the  above  equation  to  be  PqU^  —  u^  and 
also  use  the  definition  of  the  projection  Pq  ,  we  have 

B{PS^  -u^,w-Ihw)  =  {v.P^^-u^), 
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where  ///  is  the  usual  interpolation  operator  at  the  H-\eve\  nodal  points 
in  ri.  By  the  well  known  result  in  finite  element  approximation  theory,  we 
have  that 

II  w  -  Iffw  ||//^(n)  <   CH''  I  w  |//i+7(n)  <   CH'^  \\  v  Wl^^q)  . 

Therefore,  we  have 

I   B{P,^U'  -  U\  W  -  Ihw)   \<   CH-'   II   V   |U.(n)||   P,^'^  -  u'   ||^:(n)    . 

By  substituting  this  estimate  into  (2.9),  we  obtain  the  desired  proof  of  (a). 
(b)  By  the  boundedness  assumption  on  the  bilinear  form  iV(-,-),  we 
have 

I  N{PS\u'-P,^u')  \<  C\\  pS'  Wmn',)  (II  "'  L(n:)  +  II  ^.^"^  lL(n;))- 
It  is  easy  to  see  that 

II  P.^u"  ||^,n;)  =  B{pS\pS') 
=  B{u\P^u')  <  c\\u^\\^^^.^\\pS'Uuy 

which  implies  that 

II  PF^'  \W(U[)  <    C\\u'  |U(„;)  . 

The  proof  of  (b)  follows  from  combining  the  above  estimates  and  Poincare's 
inequality.  □ 

Remark:  It  is  easy  to  verify  that 

TV  N 

t=0  t=0 
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Using  Lemma  2.3,  we  can  estimate  the  right  hand  side  of  the  above  expres- 
sion by 

OT-f:i|P.u''|U,„;)||uML4,«;)- 
1  =  0 

Applying  Cauchy's  inequality  and  estimate  (1.9),  we  have  that  the  above 
expression  is  bounded  by 

1=0  '  1=0 

=  CH''B{u\Pu'') 

We  take  Ho  =  {l/Cy^\  Thus  i{  H  <  Ho,  then  6  =  CH^  <  1, 
therefore,  the  assumption  for  Theorem  1.3  part  (3)  holds. 

The  independency  of  the  mesh  parameters  can  be  accomplished  by  ex- 
amining all  the  constants  appeared  in  the  above  discussion.  The  proof  for 
Theorem  2.1  is  thus  completed. 

2.3      Iterative  substructuring  method  for  the 
problems  in  R~ 

2.3.1      An  algorithm  and  the  main  results 

1.   Basic  domain  decomposition  and  the  projections 

Recall  that  in  the  previous  section,  we  extended  each  iif-level  substruc- 
ture Cti  to  a  large  region  Q,,  on  which  we  define  our  projection  and  the 
corresponding  subproblem.  {fi,}  is  the  basic  decomposition  of  the  domain. 
Now,  instead  of  extending  each  region,  we  combine  each  pair  of  adjacent 
substructures  to  form  our  basic  decomposition  of  Ct.  We  need  to  introduce 
some  notations.  Denote  by  Tij  the  edge  which  is  common  to  two  adjacent 
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substructures  17,-  and  flj  and  let  fi.j  =  QiU  r,j  U  Q,j.  Then  {fi,j}  forms  the 
basic  decomposition  of  Q,  on  which  we  will  define  our  projections  and  this 
will  eventually  lead  to  an  iterative  substructuring  algorithm.  For  simplicity, 
we  denote  fioo  =  ^  and 

A     =  {ii,j)  1  a:i,Xj  £  A    ,x,,Xj  adjacent,  or  i  =  j  =  0}. 

Associated  with  each  Q,j,  (i,  j)  6  A^,  we  define  a  finite  element  subspace 

V,';  =  H',{n,,)nV\  for(i,j)7^(0,0). 

For  {i,j)  =  (0,0),  we  set  Fqq  =  ^^ •  ^^  i^  easy  to  verify  that 

(',J)€A£ 

We  note  that  there  is  a  big  difference  between  the  decomposition  {fi,} 
and  {^tj}-  In  the  first  case,  we  assume  that  the  overlap  is  of  order  i7, 
therefore  for  each  node  x,  G  A^,  there  is  a  ball  0(x,)  with  diameter  of 
order  H,  centered  at  x,,  such  that  0(x,)  C  ^^  for  all  ^k  with  x^  as  one  of 
its  vertex.  It  follows  from  this  property  that  we  were  able  to  use  a  local 
averaging  technique  to  obtain  to  an  optimal  estimate.  In  the  second  case, 
this  argument  can  not  be  carried  out  because  the  overlap  is  smaller,  in  the 
sense  that  no  vertex  of  A^  belongs  to  any  of  the  fi.j,  {i,j)  ^  (0,0). 

We  denote  the  projection  from  V''  to  V^^  with  respect  to  the  bilinear 
form5(•,•)byi^f,for(^,;)G  A^. 

P^  is  defined,  as  usual,  as  the  sum  of  the  projections 

(•■.J)€AS 

This  is  the  main  operator  we  shall  study  in  the  following  sections. 
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2.  The  basic  algorithm  and  tlie  main  result 

In  this  section,  we  present  our  basic  algorithm  and  the  corresponding 
convergence  rate  estimate. 

Let  u    be  the  Galerkin  approximation  of  the  solution.  By  taking 

b  =  p^u'  =    Y.  ^.>'- 

(■■,J)GA£ 

We  can  define  our  derived  equation  with  respect  to  the  bilinear  form  B{-,  •) 
and  decomposition  {Vij}  as 

P^u'   =  b. 

Next,  we  will  give  one  of  the  main  theorem  of  this  chapter,  which  esti- 
mates the  conditioning  of  the  operator  P^ .  Compared  with  the  subspaces 
used  in  the  previous  section  for  A^I,  we  use  less  overlap  here.  This  is 
reflected  in  the  poorer  bounds  on  the  operator  P^ . 

Theorem  2.2    (1)  There  exists  a  constant  Cp  >  0,  independent  of  H  and 
h,  such  that 


"y 


(2)  There  exists  a  constant  c  >  0,  independent  of  H  and  h,  such  that 

WP^u^Wa  >  c(l  +  log(i7//i))-2||u''|U,    Vt/  G  V\ 

(S)    There   exist  a  constant  Hq   >  0,   independent  of  H   and  h,   and  a 
constant  c{Hq)  >  0,  such  that,  for  H  <  Hq, 


{u\P''u^)a  >  Cp{Ho,H,h){u\u'')A,  Vu"  e  V 

where 

c^{Ho,H,h)   =  c{Ho)il+\og{H/h))-^>0. 
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A  proof  will  be  given  in  the  next  section. 

Remarks:  (a)  If  B{-,-)  is  selfadjoint,  i.e.  A{-,-)  =  B{-,-),  then  the 
operator  P^  is  symmetric  with  respect  to  A(-,  •).  Therefore,  the  generalized 
conjugate  gradient  method  can  be  used.  Theorem  2.2  part  (1)  and  (2)  give 
the  convergence  rate  estimate. 

(b)  A  remarkable  fact  about  the  algorithm  mentioned  in  (a)  is  that  if 
we  reduce  each  subproblem  P,^  to  its  Schur  complement  on  r,_,  and  then  re- 
place the  Schur  complement  by  the  square  root  of  discrete  one-dimensional 

1/2 

Laplacian,  denoted  by  /q    ■,  then  we  obtain  one  of  the  Bramble,  Pasciak  and 
Schatz's  algorithm  [4];  cf.,M.  Dryja  and  O.  Widlund  [12]. 

(c)  In  the  early  stage  of  the  development  of  domain  decomposition  meth- 
ods, a  lot  of  algorithms  were  discussed  for  two  subregion  problems,  which 
do  not  have  very  much  interest  for  practical  applications  on  multiprocessor 
machines.  Among  those  algorithms  were  the  Neumann- Dirichlet  algorithm 
[2],  modified  Schur  complement  method  [20].  However,  we  can  apply  these 
two  region  solvers  as  preconditioners  for  the  semilocal  problems  represented 
by  P,^  and  obtain  good  algorithms  for  multiprocessor  systems. 

(d)  In  general,  P^  is  nonsymmetric  but  positive  definite,  the  GMRES 
method  is  a  generally  applicable  iterative  method  for  this  class  of  prob- 
lems. Part  (3)  shows  that  the  GMRES  method  converge  in  the  (•,  •)a  inner 
product  provided  that  the  coarse  mesh  size  is  nne  enough.  The  rate  of 
convergence  depends  mildly  on  the  mesh  parameters  H  and  h. 

(e)  This  theorem  holds  only  in  2D. 

2.3.2      The  condition  number  estimates 

We  first  quote  a  lemma,  which  plays  an  important  role  in  the  traditional 
theory  for  iterative  substructuring  algorithm.     Variations  of  this  result, 
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which  dates  back  at  least  to  1966,  are  given  in  a  number  of  papers,  see  e.g. 
Bramble  [3],  Bramble,  Pasciak  and  Schatz  [4]  or  Yserentant  [29]. 

Lemma  2.4   Let  a  be  any  value  of  u^{x)  G  V^,  with  x  6  Q,,  then 


W  -  c^WUn,)  <C{1  +  \og{H/h))  I  u'  |^,(n, 


Remark:    This  result  holds  only  in  '2D.    The  1  +  \og{H/h)  factor  which 
appears  in  Theorem  2.2  is  introduced  here. 

The  next  lemma  due  to  M.  Dryja  and  O.  Widlund  [12],  provides  a 
decomposition  of  all  functions  in  V'^  into  a  sum  of  functions  defined  in 
subspaces.  This  decomposition  is  not  uniformly  bounded  as  there  is  less 
overlap,  but  the  bound  depends  only  mildly  on  the  size  parameters  h  and 
H. 

Lemma  2.5   For  any  u''  €  V",  there  exist  t/f^  G  V^'-,  {i,j)  G  A^,  such  that 


(>,j)eA£ 


4, 


and  there  exists  a  constant  C.  independent  of  u  ,  h  and  H,  such  that 
E     ll<||^,n.,)<C(l  +  log(i//h))^||u''||^ 

(»,J)6A^ 

Proof:(See  Dryja  and  Widlund,  [12]) 

Before  we  prove  theorem  2.2  part  (3),  we  need  to  formulate  a  lemma, 
which  has  already  been  proved  in  the  previous  section.  Since  the  notations 
are  a  little  bit  different,  we  rewrite  the  lemma  here  without  proof. 

Lemma  2.6    There  exists  a  constant  C  >  0,  such  that 

(a)     I  N{P,lu\u'  -  P,^y)  I  <   CH-^  II  Po^y  lUII  n'  |U,    Wu'  G  V' 
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(6)     I  N{P^^u\u'  -  P,^u')  \<   CH\\  P,^u'  |U(n.,)||  n'  Un.,).    ^u'  e  V 

/or(z,j)  7^  (0,0)  gA^. 

Remark:  By  using  the  same  argument  as  in  the  Remark  after  Lemma  2.3, 
we  can  show  that  if  the  coarse  mesh  size  is  small  enough,  we  have 

I  Niu^P^u^)  \<  6  I  B{u\P^u'')  I,     Vu''  e  V\ 

where  0  <  <5  <  1.  Therefore,  the  assumption  for  Theorem  1.3  is  established. 
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Chapter  3 

Methods  for  paraboHc 
convection-diffusion  problems 


In  this  chapter,  we  apply  the  AgA/and  /gA/to  the  parabohc  convection- 
diffusion  problems.  A  modified  AgA/is  also  introduced.  We  consider  the 
linear  systems  of  equations  that  arise  when  using  implicit  schemes  to  ap- 
proximate parabolic  problems.  If  we  consider  the  discrete  parabolic  prob- 
lem at  a  fixed  time  level,  it  is  equivalent  to  an  elliptic  problem  with  an  extra 
time  step  parameter.  The  central  mathematical  question  is  to  estimate  how 
the  convergence  rate  depends  on  the  time  step  parameter,  especially  in  the 
case  that  this  parameter  is  relatively  large,  as  well  as  the  space  mesh  pa- 
rameters. Related  works  can  be  found  in  [21],  [15]. 

The  outline  of  this  chapter  is  as  follows.  In  section  1,  we  present  the 
parabolic  problem.  In  section  2,  we  study  the  additive  Schwarz  algorithm. 
In  section  3,  we  study  a  modified  additive  Schwarz  algorithm,  which  works 
well  only  for  parabolic  problems.  Finally  in  section  4,  we  study  an  iterative 
substructuring  method. 
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3.1      A   parabolic   convection-diffusion   prob- 
lem 

We  consider  the  following  parabolic  convection-diffusion  problem:  Find 
u{x,t),  such  that 

'  du/dt  +  Lu     =    f,  in     fix[0,T], 

u    =    0  on    5Q  X  [0,r],  (3.1) 

u{x.O)     —     h   'x)     in     Q, 
where  X  is  a  strongly  ellip..:  operator  defined  in  Chapter  2.    Q  C  i?'^  is  a 
polygonal  domain  with  boundary  dCl. 

Let  us  consider  the  corresponding  variational  problem:  Find  u{x,t)  G 
ir^(Q),  t  G  [0,r],    u{x,0)  =  uo{x)  inQ,  such  that 

(-^,v)  +  B{u,v)   =   {f,v),    yvEH'oi^), 

where  the  bilinear  form  B(-,-)  and  the  linear  functional  (/,  v)  are  the  same 
as  in  Chapter  2. 

The  existence  and  uniqueness  of  the  solution  of  the  variational  parabolic 
convection-diffusion  equation  are  well  understood,  see  [19]. 

We  use  two  types  of  time  discretization,  namely,  a  backward  Euler 
scheme  and  an  implicit  Crank-Nicolson  scheme. 

Let  Atn  be  the  n"*  time  step,  M  is  the  number  of  steps  and  XInli  ^^n  = 
T.  For  the  first  scheme,  the  time  discrete  problem  is 

with  u°(a;,  t)  =  uo{x)  and  n  =  1,  •  •  • ,  M. 
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For  the  second  scheme, 

with  u°(x,  i)  =  Uo(2:)  and  n  =  1,  •  •  • ,  Af . 

To  find  the  n^^  level  solution,  both  schemes  lead  to  the  following  varia- 
tional problem:  Given  function  g-n-i  £  L^{Q).  find  w  G  Hq{Q.),  such  that 

(u;,iO   +  tB{w,v)   =   (ff„_i,r),    Vi;  6  ^^(f^).  (3.2) 

We  call  r  the  time  step  parameter.  It  is  easy  to  verify  that  for  the  backward 

Euler  scheme 

w  =     u"-u"-\ 

r  =    Ai„, 

{9n-UV)      =      t((/,v)     -     B{U'^-\V)) 

and  for  the  Crank-Nicolson  scheme 

w  =    u"  -  u"-S 

r  =    Ai„/2, 

(5n-i,W     =     r(2(/,r)    -   i?(u"-\iO). 

The  stability  of  both  schemes  is  well  understood,  see  [19].  In  this  chapter 
we  focus  on  the  study  of  fast  iterative  algorithm  for  solving  the  resulting 
large  lineax  systems  at  each  time  step.  The  variational  problem  (3.2)  is  the 
main  subject  of  our  further  study. 

To  simplify  the  notations,  we  introduce  a  bilinear  form 

Dr{w,v)   =   {w,v)   +  tB{w,v). 

In  general,  Dri:, ")  i^  not  symmetric.  For  technical  reasons  it  is  conve- 
nient to  separate  the  symmetric  and  the  skewsymmetric  part  of  this  bilinear 
form.  We  therefore  introduce  the  bilinear  forms  as 
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which  is  the  symmetric  part,  and 

which  is  the  skewsymmetric  part. 

We  assume  that  Nt{-,-)  is  bounded  in  the  following  sense:  There  exists 
a  constant  C  >  0,  such  that 

\Nr{u,v)\<CT\\uy^^a^\\v\\mn),  yu,v  eH',{n). 

Before  we  study  the  finite  element  solution  of  (3.2),  we  need  to  establish 
some  bounds  for  the  bilinear  form  Dt{w,v).  Define  the  r-norm  as  the 
square  root  of 

II  •  WriQ)  =  11  •  Wi^in)   +  "^  II  •  ll//oi(n)  • 

The  following  lemma  gives  the  boundness  and  positive  definiteness  of 
the  bilinear  form  Dt{-,-). 

Lemma  3.1    There   exist  positive   constants  Ci,   C2,   C3   and  C4,   which  are 
independent  of  t,  such  that 

(1)  I  Dr{w,v)  |<  ci  II  w  ||,,n)||  V  11,(0),    Vu;,t;  e  Hl{^). 

(2)  Dr{w,w)  >C2\\w  ||2,„,,    Vu'  €  Hlin). 

(3)  I  AAw^v)  \<  C3  II  w  11,(0)11  V  ||„n),    Vu;,i;  €  Hl{^). 

(4)  Ar{w,w)  >  C4  II  u;  ||?(o),    Vu;  e  Hl{n). 

Proof:     (1)  Because  the  boundedness  of  B(w,v),  we  have 

I   Dr{w,v)    \<\\   W    \\l2(U)\\   V    WmU)     +  Cr    II    U;   ||H>(n)||   V   ||;^.(n)    . 

The  estimate  (1)  follows  from  Cauchy's  inequality. 
The  estimate  (2)  follows  from  the  ellipticity  of  B(-,  •). 
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The  estimate  (3)  is  a  corollary  of  (1),  and  the  estimate  (4)  is  a  corollary 
of  (2).  □ 

If  we  denote  ||u||yi^  =  \jAr{u,u),  then  we  know  from  the  above  lemma 
that  the  A^-norm  is  equivalent  to  the  r-norm.  We  shall  use  ^^-norm  in  our 
discussion. 

To  obtain  a  fully  discrete  problem,  we  discrctize  equation  (3.2)  in  space 
by  using  Galerkin's  finite  element  method  in  the  subspace  V^  C  Hq{Q). 
The  finite  dimensional  approximation  of  equation  (3.2)  reads  as:  At  each 
step  n,  find  w^  G  V^,  such  that 

DAw\v')   =   (gn-uv'),    "^v'eV.  (3.3) 

This  equation  is  the  main  subject  of  this  chapter.  We  shall  use  three 
methods  to  solve  it. 

Based  on  the  H^'^'^{Cl)  regularity  assumption  on  B{u,v),  see  (2.2),  the 
equation  (3.3)  satisfies  the  following  the  regularity  result. 

Lemma  3.2   For  any  g  6  L^{fl),  the  equation 

Dr{v,w}   =  {g,v),    yv^H'oin) 

has  solution  ?r  G  H^'^^{Q)  H  Hq{Q),  and  moreover,  there  exists  a  constant 
C  >  0,  independent  ofr.w  and  g,  such  that 

I  w  |//'+^(n)<  C/T\\g\\L2^n)- 

Proof:  Consider  the  equality 

{w,w)  +  tB{w,w)  =  {g,w). 

Since  B{w,w)  >  0,  we  obtain  that 

||u'||L2(n)  <  ||^||L2(n), 
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which  imphes  that 

Ww-gWtHQ)  <  2||^||£,2(n). 
On  the  other  hand,  we  have 

B{v,iu)  =  {{g-w)/T,v),     WveH^i^). 

By  the  assumption  (2.2),  there  exists  a  constant  C  >  0  such  that 

I  w  |//i+-.(n)<  C\\{g  -  u;)/r||£,2(fi). 


3.2      Additive  Schwarz  method  for  the  prob- 
lems in  R~  and  R^ 

3.2.1      An  algorithm  and  the  main  results 

In  this  section,  we  propose  an  additive  Schwarz  algorithm  for  the  finite 
element  equation  (3.3). 

As  shown  in  Chapter  2,  the  finite  element  space  V'^  can  be  decomposed 
into  the  sum  of  a  coarse  mesh  function  space  and  a  number  of  function 
spaces  which  are  supported  only  in  subregions  0,-,  i.e. 

We  denote  by  Pf^  to  be  the  projection  from  V'^  to  V;''  with  respect  to 
the  bilinear  form  Dr{-,-)  and  let  P^^    =   P^^  H f-  P^^ 

The  additive  Schwarz  algorithm  for  equation  (3.3)  can  be  stated  as:  At 
each  step  n,  find  w^  £  V^,  such  that 
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which  is  the  derived  Galerkin  finite  element  equation  with  respect  to  the 
bilinear  form  Dr{-,-)  and  the  decomposition  {K}-  We  will  prove  that  this 
system  is  uniformly  well  conditioned.  More  precisely,  the  condition  number 
of  P^'  will  not  change  if  we  ( 1 )  refine  the  fine  mesh  size  h  to  increase  the 
accuracy;  (2)  refine  the  coarse  mesh  size  H  so  that  more  processors  can  be 
used;  (3)  increase  the  time  step  r. 

For  the  conditioning  of  the  operator  P^^ ,  we  have 

Theorem  3.1  (1)  There  exists  a  constant  Cp  >  0,  independent  of  H ,  h, 
T,  such  that 

(2)  There  exists  a  constant  c  >  0,  independent  of  H ,  h,  t,  such  that 

(S)  If  ch,t  =  nia.x{H,  H'^JH-^ /t  +  1}  is  small  enough,  i.e.  0  <  ch,t  ^ 
Co,  then  there  exists  a  constant  Cp(co)  >  0,  such  that 

{u\P''^u')Ar    >  Cp{Co){u\u')Ar,     Vu'^   G   V". 

A  proof  is  given  in  the  next  section. 

Remarks:  (a)  In  the  case  that  the  first  order  terms  in  L  vanish,  the 
problem  is  selfadjoint.  It  can  be  seen  easily  that  the  operator  P^*"  is  sym- 
metric with  respect  to  Ar-norm.  The  standard  conjugate  gradient  method 
in  .4T-norm  can  therefore  be  used.  Theorem  3.1  shows  that  this  P^''  linear 
system  is  optimal  for  the  conjugate  gradient  method  in  the  sense  that  the 
rate  of  convergence  does  not  depend  on  the  mesh  parameters  H  and  h,  nor 
on  the  time  step  size  r. 
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(b)  Since,  in  general,  (3.4)  is  a  nonsymmetric  but  positive  definite  sys- 
tem, we  use  the  GMRES  method  to  solve  it.  The  inner  product  we  use  is 

(c)  In  general,  cq  depends  on  the  coefficients  of  the  first  order  terms  in 
L,  the  ellipticity  constant  of  Dr  ,  the  bounds  on  D^  and  also  the  geometry 
of  the  domain  Vt.  We  do  not  have  a  explicit  relation  between  cq  and  the 
skewsymmetric  coefficients  in  L.  From  the  lower  bound  proof,  which  will  be 
given  in  the  next  section,  we  know  that  as  the  skewsymmetric  coefficients 
in  L  increase,  cq  decreases. 

3.2.2      Proof  of  the  theorem 

In  this  section,  we  prove  the  theorem  stated  in  the  previous  section.  Most 
of  the  techniques  are  similar  to  those  that  were  used  for  the  time  indepen- 
dent problems.  To  complete  the  proof,  we  need  only  to  show  that  all  the 
assumption  of  Theorem  1.3  hold,  and  in  addition  that  all  the  constants  that 
appear  in  the  abstract  Theorem  1.3  are  independent  of  the  mesh  parameters 
iy,  h  and  r  as  stated  in  Theorem  3.1. 

We  start  with  some  lemmas  which  contain  most  of  the  basic  results. 

The  following  lemma  is  well  known.  It  says  that  in  the  finite  element 
space  V''  the  continuous  L^  norm  is  equivalent  to  the  discrete  L^  norm.  We 
do  not  include  the  proof. 

Lemma  3.3  There  exist  two  constants  cj  >  0  and  C2  >  0,  which  defend 
only  on  the  shape  regularity  of  the  finite  element  subdivision  ofQ,  such  that 

Eere  d  is  the  dimension  of  the  space.  The  statement  is  also  true  if  we 
replace  V^  by  V^  and  h  by  H . 
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Recall  that  when  we  considered  the  .4gi\/for  the  time  independent  prob- 
lems, a  Hq{Q.)  decomposition  lemma  played  a  very  important  role  in  the 
condition  number  estimate.  Next,  we  prove  that  the  same  lemma  holds  in 
L'^{ft}  norm.  By  combining  the  estimates  in  Hq{Q.)  and  L^{Q),  we  eventu- 
ally prove  a  decomposition  lemma  in  the  ||  ■  ||  4^  norm,  which  will  be  used 
for  the  time  dependent  problems. 

Lemma  3.4  Vu''  G  V^,  there  exist  u'l  eV/",  /  =  0, 1,  •  •  ■ ,  N,  such  that 

N 

h         sr^    h 
u     =   }_^u, 

1=0 

and,  moreover,  there  exists  a  constant  C  >  0,  such  that 

N 


J2  II  ^'^  \\h(Q)^  C  II  ^'■ 


h    ||2 


li,2(n)' 


i=0 


where  C  is  independent  of  h,  H  and  u^ . 

Proof:  The  construction  of  uf,  i  =  0, 1,  •  •  •  ,  A^  is  the  same  as  in  the 
Hq{Q,)  case.  Let  ///  be  L^  projection  operator  into  the  coarse  mesh  space 
defined  in  Chapter  2.  We  take  Ug  =  Iru^-,  w^  =  u^  —  Ih^^-,  aJid  then  set 
u,^  =  Ik{Bi{w^)),  where  the  partition  of  unity  {^,}  is  the  same  as  in  Chapter 
2. 

Because  Ihu^  is  the  X^  projection,  we  have 

II  Ihu^  \\Ihu)<\\  ^'  WIhu)  ■  (3.5) 

Next,  we  prove  the  X'^(f2)  boundness  for  the  other  uj".  By  using  Lemma  3.3 
for  each  Q-, 
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Since  |  9,  \<  1,   i  =  1,  •  •  •  ,iV,  and  {Vi\}  is  a  finite  covering  of  Q,  the  right 
hand  side  can  be  bounded  by 

Using  Lemma  3.3  again,  we  obtain 

E  II  ^-  WUu  ><  C  II  ro'  \\h(^)<  C  II  u^  lli.(n)  .  (3.6) 

1=1 

Combining  the  estimates  (3.5)  and  fT  G)  proves  this  lemma. □ 

The  decomposition  ler:  .  for  tht  |j  •  ||T(n)  norm  follows  immediately 
from  Lemma  3.4  and  the  Hl{Q.)  decomposition  Lemma  2.2.  We  simply 
state  it  as  follows: 

Lemma  3.5  \/u^  e  V^.  there  exist  uj  G  V;'',  t  =  0, 1,  •  •  •  ,  A^,  such  that 

N 
1=0 

and,  moreover,  there  exists  a  constant  C  >  0,  such  that 

f:ii"Mii(.',<cii.''iii, 


r(n    ) 

1=0 


where  C  is  independent  of  h,  H  and  u 


h 


Lemma  3.6  For  any  w  E  H^ifl)nH^-^'^{Q,),  there  exists  a  w"  E  V" ,  such 
that 

\\w  -  w"\\r(u)  <  CHWH^  +  T  I  w  |H'+.(n), 

where  C  is  independent  of  t ,  w  and  H. 
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Proof:  For  any  given  w,  let  w^  be  the  solution  of 

a(w^,v)   =   a(w,v),     \/v  E  V" , 

where  a{u,v)  =  /q  VuVvdQ..  By  the  classical  finite  element  approximation 
theory,  we  know  that 

\\w  -  w^\\hi(q)  <  CH''  I  w  |//i+^(n), 

\\w  -  w"\\l2(q)  <  CH^'^''  I  w  |Hi+7(n)  . 
Hence,  we  have 


Ik  -  w"\\r^Cl)   <  CH-'VH^  +  T   I   W   |;/i  +  .(n)   • 

D 

Recall  that  Pq  ''  is  the  projection  from  V'^  into  V^  with  respect  to  the 
.Dr{-,  ■)  norm,  thus  a  H^  projection  in  some  sense.  We  can  therefore  expect 
the  L^  approximation  to  contribute  an  extra  factor  H.  The  L^{^)  estimate 
is  given  in  the  following  lemma. 

Lemma  3.7    There  exists  a  constant  C  >  0,  which  is  independent  of  both 
H  and  t,  such  that 

lit/'' - p^^u'^Wm^)  <  CH'^{y/ipT^)/T\\u%(^),  \/u'  E  y\ 

Proof:  We  first  establish  a  bound  for  the  r-norm.  It  is  easy  to  verify  that 

Using  Lemma  3.1  for  both  sides,  we  get  the  estimate 

||u''-Po^^u''||.(0)<Cl|u''||,(n). 
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The  ^o(^)  (or  r-norm)  estimate  and  the  Nitsche  trick  give  the  L^(f)) 
estimate.  We  have 

\\u'  -  Po'^'u'U.^n)  =       sup       '- —^ ^. 

For  any  fixed  v  G  L^(fi),  we  form  an  auxihary  variational  problem:   Find 
w  e  ifi+^(n)  n  H^iQ),  such  that 

By  Lemma  3.2,  we  have  |  '.     -/i+-'(n,l'   -"'/■''Ili'llL^in). 
We  take  cp  =  u^  -  P^'  l      ^  H^(n)\  then 

{u^  -  P^'u\v)  =  Dr{u^  -  P^'u^,w). 

Take  w^  to  be  the  V^  approximation  of  w  obtained  in  Lemma  3.6.  Since 
w^  €  V^ ,  we  have 

\Dr{u''  -P^^U^,W)\       =       \Dr{u''  -P^^U\W-W")\ 

<  C\\u^-P^^u%(u)\\w-w^\Un^ 

<  CH■^^/JPT^\\U%(^^  I  W  |;,,+,(n)  . 

Combining  the  above  results,  we  have 

llu''  -  P^^u'^Wm^^  <  CH-^{VlPT^)/T\\n%(u). 

which  is  the  desired  result.  □ 

In  the  next  lemma,  we  estimate  the  energy  contributed  by  the  skewsym- 
metric  part  Nr{--,  •)•  We  show  that  it  is  in  fact  a  lower  order  term  compared 
with  the  symmetric  part,  and  hence  can  be  controlled  if  the  coarse  mesh 
size  is  fine  enough. 
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Lemma  3.8  If  xnax{H ,  H''  J H'^ / t  +  1}  is  small  enough,  there  exists  a  con- 
stant 0  <  6  <  1,  such  that 

Proof:  Since  Nr[PF'u^,PF'u^)  =  0  ,  it  is  easy  to  verify  that 

N 
1=0 

Therefore,  we  need  only  to  estimate  the  right  hand  side  of  the  above  in- 
equahty.  Since  the  coarse  mesh  projection  P/'  is  special,  we  will  consider 
it  separately  using  the  result  of  Lemma  3.7. 
(1)  If  i  =  0,  by  Lemma  3.7 

\N,{PF^u\u^-PF^u')\     <     CT\\PF^u'\\H.^^)\\u'^-Po''^u'^\\mu) 

<    CH-VWTT\\P,''^u'^\\Hi^ny\\u%^n). 

It  is  easy  to  see  that 

\/r|K'i|//.(Q)  <  ||r'|,-(n). 

Hence, 

I   NriPo'^'u'^u'   -  P^^U')    |<   CH-^m/T+l\\PF^U%(^^\\u%(uy 

From  the  definition  of  Pq'',  we  have 

If  we  apply  Lemma  3.1  for  both  sides,  we  obtain 
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Therefore,  the  first  term  can  be  bounded  as  follows 


(2)If  i  ^  0, 

<    c^\\pp^u'\y^^.^{\\u%^^.^  +  \\PF'u%^^.^). 

The  first  factor  can  be  estimated  by  using  Poincare's  inequality.     Since 
P,  ''u^  G  Hl{Q.\),  the  diam-'^er  of  fi,       of  order  if,  and  we  have 

ll^.''^"''iiL^,Q:)<ci;ri|p.^^u''ii^.(„;). 

By  the  definition  of  P[^% 

Using  Lemma  3.1  for  both  sides,  we  obtain 

Combining  these  inequalities,  we  obtain 

I  Nr{PF'u\u'  -  PP^u')  \<  CH\\uX^^,y 
Putting  the  results  in  (1)  and  (2)  together,  we  have 

I  NriP'^^u^  u'')\<C  max{if,  ^Vi^Vr  +  1}  ^  hX^n'r  (3.7) 


1=0 


which  also  holds  in  ^Ir-norm.  It  is  ea^y  to  see  that 

N 


j:\\pF'A\i^o'j  =  D.{u\p^^u^). 


t=0 
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By  using  the  decomposition  lemma,  we  have 


<  crLo\\pru'^\Uin:)\K\\A.in:) 


<  c^z^Au^r,^,,.,^zLo\lPF^uX^,,y 


Therefore,  we  obtain 

1=0  ' 

Hence 

\\u'\\\^<CDAu\P''-x'').  (3.8) 

We  complete  the  proof  by  combining  estimates  (3.7)  and  (3.S).  □ 

The  proof  of  Theorem  3.1  can  be  accomplished  by  using  Theorem  1.3 
ajid  the  lemmas  in  this  section. 

3.3      Modified  additive  Schwarz   method  for 
the  problems  in  R~ 

3.3.1      An  algorithm  and  the  main  results 

In  this  section,  we  propose  a  modified  version  of  .4s.^/by  dropping  the  coarse 
mesh  space  Vq,  which  represents  the  global  information  transportation.  We 
show  that  in  some  situations  the  global  function  space  is  not  necessary  for 
fast  convergence.  This  is  true  only  for  parabolic  equations. 

Let  us  define  P^^  =  P/^'  +  •  •  •  +  P;^^  where  the  P,^'  are  the  same  as 
in  the  previous  section.  Note  that  we  excluded  P^''. 

We  use  the  seime  algorithm,  presented  in  the  previous  section,  with  P^*" 
replaced  by  P^''.  About  the  conditioning  of  the  operator  P^"",  we  have  the 
following  theorem. 
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Theorem  3.2   (1)  There  exists  a  constant  Cf,  >  0,  such  that 

\\P'''u'\\Ar<Cp\\u'\U.,    Vu''ei^\ 

(2)  There  exists  a  constant  c  >  0,  such  that 

WP'^-u'iUr  >  c(i  +  r/H'r'WAUr,   Vu''  e  v. 

(3)  If  Cffr  =  H{\  +  T / H'^)  IS  small  enough,  i.e.  0  <  c//,.  <  cq,  then 
there  exists  a  constant  Cp(co)  >  0.  such  that 

A  proof  is  given  in  the  next  section. 

Remarks:  (a)  For  symmetric  problems,  part  (1)  and  (2)  of  this  theorem 
shows  that  if  the  factor  t/H^  is  small,  the  dropping  of  the  coarse  mesh 
space  does  not  lead  to  slow  convergence.  This  suggests  we  might  use  the 
modified  -45A/in  the  case  that  we  have  a  relatively  small  time  step  or  large 
substructures. 

For  nonsymmetric  problems,  in  order  to  obtain  the  fast  convergence,  we 
need  both  H  and  r/H'^  to  be  small.  This  implies  that  we  cannot  choose  r 
and  H  independently. 

(b)  This  theorem  is  true  only  in  R^.  In  higher  dimension,  the  fine  mesh 
size  h  enters  our  bounds. 

3.3.2      Proof  of  the  theorem 

The  proof  of  Theorem  3.2  can  be  accompHshed  by  using  the  results  in 
Theorem  1.3,  a  new  decomposition  lemma  given  below  and  the  new  estimate 
of  the  skewsymmetric  part.  We  begin  with  the  decomposition  lemma. 
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Lemma  3.9  Vu''  G  V'',  there  exist  u'^  G  V^^,  z  =  1,  ■  •  •  ,  A^,  such  that 

N 


n'   =   E"f 


1=1 
and,  moreover,  there  exists  a  constant  C  >  0,  such  that 

N 


i  =  l 

where  C  is  independent  of  h,  H  and  t. 

Proof:  We  first  construct  the  decomposition,  then  we  do  the  estimates  in 
both  Hl{Q.)  and  L'^{Q.)  norm. 

Let  {6,,i  =  1,---,A^}  be  the  partition  of  unity  defined  in  Chapter  2. 
Denote 

For  each  substructure  fi,,  we  have 

where  the  sum  is  taken  over  all  adjacent  pairs  of  nodal  points  x/  and  x^  in 
n,.  Let  A'  C  ^i  be  a  single  element  and  x/,  x^^  G  A'.  Denote 

We  then  have 

(^.u'')(x,)-(^,u'')(x^)   = 

which  can  be  bounded  from  above  by 

C{h/H max{\  u\x)  \}+  \  u''(x/)  -  u\x,^)  |). 

K 
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By  squaring  this  estimate,  using  the  triangle  inequaHty  and  summing 
over  all  K  C  i^[,  we  obtain 

Using  Poincare's  inequality,  we  have 
Therefore 


Hence,  we  obtain 


<       (^17-2    I       h   12 


1=1 


Combined  with  the  £^  result  given  in  Lemma  3.4  and  the  definition  of 
the  r-norm,  we  complete  the  proof.  □ 

In  the  next  lemma,  we  estimate  the  energy  contributed  by  the  skewsym- 
metric  part  Nt{-,-). 

Lemma  3.10   There  exists  a  constant  C  >  0,  which  is  independent  of  h,H 
and  T,  such  that 

I  Nr(u\P^^u'')  \<  CH{1  +  T/H^)Dr{u\p'^^u''),    Vu'*  G  V\ 

The  proof  can  be  obtained  by  using  the  proof  of  Lemma  3.8  part  (2) 
and  replacing  the  decomposition  Lemma  3.5  by  Lemma  3.9. 
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3.4      Iterative  substructuring  method  for  the 
problems  in  R^ 

3.4.1      An  algorithm  and  the  main  results 

1.   Basic  domain  decomposition  and  projections 

We  assume  that  the  decomposition  of  Q  is  {Q.^,  (z,j)  £  A^},  which 
is  the  same  as  the  one  introduced  in  Chapter  2.  The  corresponding  V'^ 
function  space  decomposition  is  also  the  same 

v  =    E  ^  • 

(■.j)€A£ 

Denote  by  P,^'"  to  be  the  projection  from  V^  to  subspace  V/-  with  respect 
to  the  bihnear  form  Dr{-,  ■)■,  which  vas  defined  in  Chapter  2. 
The  mapping  P^^  :    V"  — y  V^  is  defined  as 

P^^  =     E    ^.?^ 

(■.j)ga£ 

which  is  one  of  the  main  operators  we  shall  deal  with  in  this  chapter. 
2.  The  basic  algorithm  and  the  main  results 

In  this  section,  we  first  introduce  our  iterative  substructuring  algorithm 
for  the  parabolic  convection-diffusion  equation.  Then,  we  state  the  theorem 
concerning  the  rate  of  convergence  of  the  algorithm.  The  proof  of  the 
theorem  will  be  given  in  the  next  section. 

Let  us  denote 


^Dr       ^      pDr^=  J2         P!:U. 


(i,j)6A£ 


Our  basic  iterative  substructuring  algorithm  for  solving  equation  3.3  is  as 
follows: 
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Basic  algorithm:  Solve  the  following  linear  system  of  equations 

by  the  GMRES  method  using  the  inner  product  [•,]  =  (■,-)Ar- 

We  next  present  some  hounds  concerning  the  conditioning  of  the  opera- 
tor P^"".  The  rate  of  conveigrnce  of  the  above  algorithm  can  be  estimated 
by  using  these  bounds.  1        is  denote 

Theorem  3.3   (1)  There        -is  a  constant  Cp  >  0,  independent  of  H ,  h, 
T,  such  that 

(2)  There  exists  a  constant  c  >  0,  independent  of  H ,  h,  r,  such  that 


,rh 


(3)  If  CH,T,h    =   ma.x{H,H''JH'^/T  +  l}u   is  small  enough,    i.e.     0    < 
c//,T,/i  ^  Co,  then  there  exists  a  constant  Cp(co)  >  0,  such  that 

{u\P^^u'')a.    >CpU-\h,H.T)iu\u')Ar,     Vu''   G   V". 

3.4.2      Proof  of  the  theorem 

In  this  section,  we  prove  the  theorem  stated  in  the  last  section.  Again, 
we  need  to  prove  that  the  assumptions  made  for  the  Theorem  1.3  hold 
and  also  show  how  the  constants  depend  on  the  mesh  parameters  h,H 
and  r.  The  proofs  look  more  complicated  because  of  the  lack  of  overlap 
between  the  substructures.  But  the  idea  behind  the  proof  is  similar.  We 
first  establish  a  bounded  decomposition  lemma  for  V'^  functions.  Note  that 
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the  decomposition  is  no  longer  uniform  in  ||  •  ||o,  norm;  the  bounds  depend 
mildly  on  the  size  parameters  h^H. 

We  rewrite  Lemma  2.4  in  the  way  that  we  shall  use  later  in  this  section. 

Lemma  3.11   Let  Q,  be  a  substructure,  u^  r  V','',  denote 

area{fl,)  Jci, 
Then 

Lemma  3.12    For  any  u''  G  V^,  let 

h  V^  h 

be  the  decomposition  constructed  in  Lemma  2.5.    We  then  have 

E  ll"iHL(f7„,<c(i|u''iii.(.,  +  iJ^(i  +  iogi?A)llu''ii^,(^,), 

where  the  constant  C  >  0  is  independent  of  H  and  h. 

Proof:  For  a  given  u^  e  V",  we  begin  with  the  estimate  of  the  L^(n) 
norm  of  Uq^  -■  Ih^'^-  Let  us  consider  one  substructure  at  a  time.  Assume 
that  Cli  has  vertices  Ti,T2,r3  and  denote 


Q,    =    -r—  /     u'dQ. 

area! 


a(fi,)  ^n, 
Then,  we  have 
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Since  the  function  ///u''  —  a,  is  linear  in  the  region  f2,,  a  straightforward 
calculation  shows  that  the  L^  norm  in  i7,  can  be  bounded  by 

Using  Lemma  3.11,  the  above  sum  is  bounded  by 

CF^(l+logF/M|uM^>(f,,). 
To  bound  the  ||Q,||z,2(n,)  term,  we  need  to  estimate  the  following  integral 

Applying  Holder's  inequality,  this  is  bounded  by 

area(a)( ^)'  /    1^1^  •  /  (u')'dn    =    \\u'\\l,^^^y 

Therefore,  we  obtain 

ll^/z^'llLiQ.)  <  Cdlu^lli.jo,)  +  H\l  +  log  H/h)  I  u'  |^,(„,)). 

By  summing  this  inequalities  for  i  —  1,  ■••,.¥,  and  using  Poincare's  in- 
equality and  replacing  the  H^  seminorm  by  the  H^  norm,  we  obtain  the 
estimate  for  Uqq  as, 

ll^i/t^'lli^ff.)  <c(||u''|ii.,„,  +  if^(i  +  iog^/h)llu'>ll^,,„)). 

Now  consider  the  case  (z,i)  ^  (0,0).  Recall  the  construction, 


Let  us  denote 


Uf.      =     Iki6:Au'   -  Ihu"")). 


Aj  =  a,nA\ 
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By  Lemma  3.3,  we  have 

Since  I  dij  |<  1,  and  x^-  is  nodal  point,  the  interpolation  operator  Ih  can  be 
removed.  Therefore  the  right  hand  side  can  be  bounded  by 

Ch'    Yl    (iu'-lHu')(x,))\ 

If  we  sum  over  all  (?,  j)  G  A^\(0,  0),  taking  into  account  of  the  fact  that  at 
each  nodal  point  the  function  value  contribute.;  at  most  three  times,  then 

E    ii"inL(C7.,)  <  ch'  y:  (C'^'  -  iHn'){x,))\ 

(ij)eA-E\(0,0)  x/iGA'' 


Using  Lemma  3.3  again,  we  can  bound  the  right  hand  side  by 


C\\u^   -  Ihu''^^^ 


which  can  be  bounded  by 

CWn'Wh^n)  +  H\l+ log  H/h)\\uX.^^y 
The  proof  is  completed  by  adding  overall  {i,j).  Q 
Lemma  3.13   For  any  u^  £  V^,  let 

-'  =    E  4 

('.j)eA£ 
be  the  function  decomposition  constructed  in  Lemma  2.5,  then  there  exists 
a  constant  C  >  0,  independent  of  H,  h,T,  such  that 

E     H\U,,)^Cu{h,H,T)\\urriay 
(',j)eA£ 
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Proof:     The  proof  is  accomplished  by  combining  the  Lemma  3.12  and  the 
Lemma  2.5.  Since 

E(.„)eA-||uf,||?(n,,     <     C{\\u'\\h^^^  +  HHl  +  \ogH/h)\\u'^\\l,^^^ 

+Til  +  \ogH/h)\\u'^\\l^^^^). 
which  is  bounded  bounded  from  above  by 

o 

Lemma  3.14    There  exists  a  constani  C  >  0.  independent  of  h,H  and  r, 
such  that 


I  Nr{u\P'^^u^)  |<  C msix{H,  H''  y/m/r  +  l}u{h,  H,T)Dr{u\  P^^u^),    \/u^  6  V 

Proof:  It  is  easy  to  see  that  the  coarse  mesh  projection  Pqo''  is  identical 
to  the  coarse  mesh  projection  Pq  "^  defined  in  section  2  of  this  chapter.  It 
follows  from  the  proof  of  Lemma  3.S  part  (1)  and  (2),  that  we  have 

I  NriP'^'u^u')  \<  Cmax{H,H-'^HyT  +  1}     ^     \\n%n„y      (3-9) 

By  using  the  same  argument  as  when  we  proved  equation  (3.8)  and  the 
decomposition  lemma  3.13,  we  have 

ll^'lli    <  CuiK  if,  T)Driu\  P^'^U').  (3.10) 

D 

The  proof  of  Theorem  3.3  is  completed  by  using  Theorem  1.3,  the  lem- 
mas in  this  section  and  the  assumption  that  Tna.x{H,  H'*JH'^/t  +  l}u{h,  if,  r) 
is  small  enough. 
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Chapter  4 
Numerical  results 

4.1      Stationary  convection- diffusion  problems 
4.1.1      The  model  problems 

In  this  section,  we  present  some  model  problems  that  we  use  to  test  the 
additive  Schwarz  algorithm  and  the  iterative  substructuring  algorithm  for 
the  stationary  convection-diffusion  equation.     We  consider  the  following 
linear  second  order  elliptic  problem  defined  Oii  Q  =  [0, 1]  x  [0, 1]  C  R^, 
d  ,^du^       d  ,    du  .         du       ^du  . 

with  the  homogenous  Dirichlet  boundary  condition.  The  coefficients  are 
specified  as  follows. 

Example  0.  E,  —  \-.  rj  =  1  and  a  =  /?  =  7  =  0.  This  is  a  selfadjoint 
problem.  We  shall  use  to  test  the  iterative  substructuring  algorithm.  /  is 
chosen  so  that  the  solution  has  the  form  u   =   xe'^^sin{Trx)sin{TTy). 

Example  1.  ^  =  I  +  x^  +  y^,  i]  =  e^^  a  =  5{x  +  y),  ^  =  1/(1  +  x  +  y) 
and  7  =  0.  u  is  the  same  as  in  Example  0. 

Example  2.  The  coefficients  are  chosen  as  ^  =  cr,  77  =  cr,  a  =  1,  /3  =  1 
and  7  =  1.  <7  will  be  specified  later,  u  is  the  saxne  as  in  Example  0. 
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The  ^-level  subdivision  of  Q  is  introduced  in  a  simple  way  shown  in 
the  following  picture.  We  further  divide  each  subregion  in  the  same  fashion 
into  /i-level  triangles,  which  are  not  shown  in  the  picture. 

Y 


Figure  1.  H— level  subdivision 

Let  us  introduce  some  notations,  ite  denotes  the  total  number  of  itera- 
tions required  for  the  GMRES  algorithm,  max  err  denotes  the  maximum 
error  between  the  numerical  solution  of  the  Galerkin's  equation  and  the 
exact  solution  of  the  continuous  equation,  ovip  measures  the  size  of  the 
overlap,  that  is  a  integer  multiple  of  the  fine  mesh  size  h.  In  our  programs, 
all  the  subproblems  are  solved  exactly  by  the  band  solver  from  LINPACK. 
The  stopping  criterion  for  the  GMRES  method  is  ||  r,  ||^  /  ||  tq  ||<  10~^, 
where  r,  is  the  i"'  step  residual  defined  in  Chapter  1.  The  programs  were 
run  in  single  precision  on  a  CONVEX  C-1  computer  at  NYU. 

4.1.2      Tests  of  the  additive  Schwarz  method 

Based  on  the  2-level  subdivision  of  f2,  we  introduce  our  basic  subregion  Cl- 
obtained  by  extending  each  triangle  f2,  to  a  larger  triangle  such  that  each 
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side  is  parallel  to  its  corresponding  side,  and  does  not  cut  through  any  h- 
level  triangles.  We  assume  the  extension  of  each  side  of  the  i7,  is  uniform 
in  the  sense  that  they  all  extend  the  same  number  of  /i-level  triangles.  See 
the  following  picture.  We  cut  off  anything  which  is  outside  f2. 

Y 


Figure  2.  extended  subregions 

We  note  that  some  subregions  Q^  that  intersect  with  the  boundary  dO. 
are  no  longer  triangles.  A  lot  of  regions  v. i;.h  different  shapes  are  intro- 
duced even  if  the  original  i7,  has  a  nice  shape.  This  causes  difficulties  with 
applying  fast  solvers  in  practical  applications.  There  are  some  ways  around 
it,  but  more  on  the  theory  and  the  numerical  experiments  still  need  to  be 
carried  out.  For  example,  we  can  replace  the  Q-s  by  rectangles  with  ap- 
proximately the  same  size,  so  that  the  FFT  can  be  used  as  a  subproblem 
solver. 

Example  1.  In  this  example,  we  plan  to  test  the  effectiveness  of  the 
algorithm  for  different  mesh  sizes  h  and  H  with  a  fixed  overlapping  factor 
a.  But,  since  the  overlap  must  set  to  be  a  integer  multiple  of  the  fine  mesh 
size  h,  it  turns  out  that  we  cannot  always  use  the  same  overlapping  factor 
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for  all  h  and  H.   In  the  following  table,  \{  H  =  1/3,  we  set  the  overlap  to 
be  approximately  OAH ,  and  if  if  =  1/5,  to  be  0.33H . 


h 

H 

ite 

H 

ite 

max  err 

1/15 

1/3 

14 

1/5 

12 

5.0  X  10-3 

1/30 

14 

1/5 

14 

1.3  X  10-3 

1/45 

1'^' 

15 

1/5 

15 

5.8  X  10--* 

1/60 

15 

1/5 

15 

3.0  X  10-^ 

For  fixed  H.  the  opti  y  in  r.   :an  be  seen  clearly,  i.e.    ite  changes 

very  slowly  if  we  use  the  finer  mesh  size  h.  However  if  we  change  H,  the 
number  of  iteration  required  to  achieve  the  same  accuracy  does  change. 
That  is  related  to  the  parameter  6  in  Theorem  1.3  which  is  increased  if  the 
coarse  mesh  is  further  refined  subject  to  the  constraint  H  <  Hq.  This  can 
be  seen  more  clearly  for  problems  with  higher  Reynold's  number  as  in  the 
following  examples. 

In  the  theory  for  additive  Schwarz,  we  assume  that  the  distance  between 
the  interior  boundary  of  Q,  ana  the  boundary  of  J7,  is  to  be  bounded  from 
below  by  a  constant  a  times  the  coarse  mesh  size  H.  Now,  in  the  following 
table  we  show  how  this  constant  affects  the  number  of  iterations.  We  use 
the  same  model  problem  as  above.  Set  H  =  1/5  and  h  =  1/60.  We  run  the 
same  program  for  different  sizes  of  overlap,  which  range  from  one  fine  grid 
size  to  six  grid  size.  The  results  are  shown  in  the  following  table. 


ovip 

h 

2h 

Zh 

4/i 

5h 

6h 

ite 

IS 

IS 

16 

15 

14 

13 
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Indeed,  increasing  the  overlap  reduces  the  number  of  iterations.  But  we 
have  to  note  that  this  also  increases  the  size  of  the  subproblems  and  the 
cost  per  iteration  is  increased  also.  Further  studies  on  how  to  determine 
the  optimal  overlap  so  that  the  parallel  CPU  time  or  the  serial  CPU  time 
can  be  minimized  are  definitely  needed.  Some  discussions  about  this  issue 
for  symmetric  problems  can  be  found  in  [16]. 

Example  2.  (1)  a  =  ^2/30. 

We  test  the  algorithm  with  a  higher  Reynold's  number.  We  first  take 
h  =  1/45  and  vary  the  parameter  H. .  It  can  ^  seen  from  the  following  table 
that  when  H  <  Hq  =  1/9  the  number  of  i:  aiions  does  not  depend  on  H 
any  more,  therefore,  the  algorithm  is  seen  lo  oe  optimal.  In  the  following 
tables  the  overlapping  factor  is  between  1/3  and  1/4, 


H 

1/3 

1/5 

1/9 

1/15 

ite 

+50 

37 

20 

11 

ovlp 

5h 

3/7 

2h 

Ih 

As  pointed  in  Theorem  2.1,  the  observed  parameter  Hq  should  not 
change  if  we  change  the  mesh  size  parameters.  This  is  demonstrated  in 
the  following  table. 


H=l/10 

h=l/S0,ovlp=2/7 

;?-l/100,ovlp=3/i 

/i=l/120,ovlp=4/i 

ite 

19 

20 

20 

(2)  a  =  \^/100 

In  this  example,  we  increase  the  Reynold  number  to  100.0.  It  can  be 
seen  that  by  choosing  a  suitable  H  <  Hq,  -we  can  still  control  the  total 
number  of  iterations  taken  to  achieve  a  specified  accuracy.  By  runing  the 
same  example  on  a  coarser  fine  mesh  space  an  approximate  Hq  can  be 
obtained.    We  refine  the  mesh  so  that  the  required  accuracy  is  achieved. 
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According  to  Theorem  2.1  this  further  refinement  of  the  fine  mesh  should 
not  increase  the  number  of  iterations.  A  further  refinement  of  the  coarse 
mesh  reduces  the  total  number  of  iterations.  See  the  following  tables. 


/i=l/45,ovlp=l/i 

/i=l/60.ovlp=l/i 

/i=l/75,ovlp=2/i 

i?=l/15 

23 

23 

23 

/i=l/SO,ovlr     Ih 

/)  =  l/100,ovlp=2/i 

/i=l/120.ovlp=2/! 

if=l/20 

15 

21 

21 

4.1.3      Tests  of  the      era  substructuring  method 

Since  few  numerical  exper....ents  of  ig^l/for  symmetric  problems  have  been 
done  before,  we  first  present  some  results  for  the  symmetric  problems.  We 
then  concentrate  on  the  convection-diffusion  problems.  Compared  with  the 
A3M,  the  JgMis  easier  to  implement  because  we  do  not  need  to  determine 
the  overlap  size. 

We  take  Q,  to  be  the  same  as  for  AgA/,  and  let  Qij  be  the  union  the 
adjacent  pair  fi,  and  flj.  The  subproblems  are  solved  exactly. 

Example  0.  Since  the  operator  P^  in  this  example  is  symmetric  with 
respect  to  the  (•,-)^,  see  [11],  the  GMRES  method  is  equivcdent  to  the 
conjugate  gradient  method. 

Let  us  first  fix  H  and  vary  /i. 


h 

1/12 

1/24 

1/36 

1/48 

1/60 

1/72 

H=l/3 

12 

13 

14 

14 

14 

14 

In  the  next  set  of  examples,  we  fix  h  and  vary  if. 


H 

1/2 

1/3 

1/4 

1/5 

1/6 

1/10 

1/12 

1/15 

1/20 

h=l/60 

11 

14 

15 

14 

14 

13 

12 

11 

10 
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As  predicted  by  the  theory,  the  number  of  iteration  depends  only  mildly 
on  the  mesh  parameters  h  and  H.  The  restriction  on  the  coarse  mesh  size 
does  not  apply  in  the  symmetric  case.  We  will  see  later  that  refinement  of 
the  coaxse  mesh  can  considerably  reduce  the  number  of  iterations  for  non- 
symmetric  problems,  however,  this  is  not  true  for  the  symmetric  problems. 
As  an  example,  we  repeat  the  first  set  of  test'^  shown  above  with  a  finer 
coarse  mesh, 


h 

1/12 

1/24 

1/36  1  "  48 

:/60 

1/72 

H=l/4: 

12 

14 

14  1    - 

15 

15 

Example  1.  This  is  a  nonsymmetric  •  blem.  However,  the  nonsym- 
metric  part  is  not  large  compared  with  the  symmetric  part.  In  the  following 
table  we  show  how  the  number  of  iterations  depend  on  the  mesh  parame- 
ters. 


h 

1/15 

1/30 

1/45 

1/60 

1/75 

1/90 

1/105 

H=l/3 

16 

IS 

19 

20 

20 

20 

21 

H=l/5 

14 

17 

IS 

18 

IS 

18 

18 

We  can  see  clearly  that  the  refinement  of  ^he  coarse  mesh  can  reduce  the 
total  number  of  iterations,  which  is  not  true  for  the  symmetric  problems; 
cf.  Theorem  2.2. 

Example  2.  (1)  a  =  n/2/30. 

For  the  problems  with  higher  Reynold's  number,  the  effect  of  the  coaxse 
mesh  size  becomes  clearer.  We  first  use  h  =  1/45  to  run  a  series  tests  for 
different  H.  Results  are  shown  in  the  following  table, 


H 

1/5 

1/9 

1/15 

/i=l/45 

41 

29 

18 

which  tells  us  that  the  proper  coarse  size  is  about  H  =   1/15.     This  is 
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confirmed  by  the  following  tests  for  different  h. 


H 

1/6 

1/10 

1/12 

1/15 

1/20 

h=l/60 

37 

23 

20 

20 

13 

/i=l/120 

43 

28 

24 

20 

17 

(2)  a  =  v^/100. 

In  order  to  maintain  atively  small  number  of  iterations  for  problems 

with  high  Reynold's  nun  we  need  to  select  a  finer  coarse  mesh,  which 

depends  strongly  on  the  Re  lold's  number  but  not  the  fine  mesh  size.  This 
fact  suggests  we  should  mine      .e  proper  coarse  mesh  by  running  a 

series  of  tests  for  problem.  .  smaller  size  with  the  same  Reynold's  number. 
Once  the  coarse  mesh  size  is  determined,  we  can  further  refine  our  fine  mesh 
so  that  the  desired  accuracy  can  be  obtained.  The  number  of  iterations 
grows  mildly  as  h  decreases. 


h 

1/30 

1/45 

1/60 

1/75 

1/90 

H=l/15 

20 

24 

30 

34 

36 

h 

1/6C   .  SO 

1/100 

1/120 

1/140 

H=l/20 

20    23 

26 

28 

30 

4.2      Parabolic  convection-diffusion  problems 

In  this  section,  we  consider  the  parabolic  convection-diffusion  equations 
defined  on  the  unit  square  [0, 1]  x  [0, 1]  C  R^,  which  have  been  discretized 
by  either  the  backward  Euler  scheme  or  the  implicit  Crank- Nicolson  scheme 
in  time  and  by  a  Galerkin  finite  elenient  method  in  space.  We  present 
numerical  results  for  computing  the  solution  at  a  fixed  time  level.  In  the 
following  examples,  the  elliptic  part  of  the  parabolic  equation  are  chosen 
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to  be  the  same  as  in  section  1  and  all  the  space  discretizations  that  we  use 
in  this  section  are  the  same  as  those  introduced  in  the  last  section. 
We  assume  that  the  time  step  has  the  form 

where  e  >  0  is  called  the  time  step  tolerance.  The  main  issue  here  is  to 
determine  how  the  rate  of  convergence  (or  r.umber  of  iterations)  depends 
on  e.  In  the  case  e  =  2.0,  the  stiffness  ma-  x  is  well  conditioned  in  the 
sense  that  the  condition  number  does  not  d'  .d  on  the  mesh  parameters  h 
and  T.  In  general,  the  smaller  e  is,  i.e.  the  i  tr  the  time  step  is,  the  more 
ill-conditioned  is  the  problem.  If  e  is  too  sniai..  ihe  accuracy  of  the  solution 
is  completely  lost.  We  do  include  some  test  results  using  unusually  large 
time  step  r,  compared  with  the  space  mesh  size  h,  to  show  the  robustness 
of  the  algorithm.  In  our  experiments,  we  assume  that  0.25  <  e  <  1.5. 
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4.2.1      Tests  of  the  additive  Schwarz  method 


e 

1.5 

1.0 

0.5 

0.25 

Example  0 

h=l/45.H=l/3,ov\p=5h 

10 

12 

12 

13 

h=l/60.i/  =  l/3,ovlp=6/i 

10 

12 

13 

13 

/2  =  l/60.i/  =  l/10.ovlp=2/! 

10 

11 

12 

12 

/^=l/4^i7  =  l/3.ovlp=/^ 

15 

IS 

IS 

19 

Example  1 

h  =  l/:     H=l/3,ovlp=2h 

12 

13 

14 

14 

h=l/       H=l/3.ov\p=4:h 

12 

14 

14 

14 

h=l/-     ¥=l/3,ovlp=6/7 

12 

14 

15 

15 

h^l/C    H=l'''    -lp=Sh 

11 

14 

15 

15 

h  =  \l       ''-1           \p=\h 

11 

12 

12 

12 

h-\,        .=l/u..    \p=2h 

12 

14 

14 

14 

h=\l^:..H  =  l/h.o\\p=3h 

12 

13 

14 

15 

/i=l/60.i7=l/5,ovlp=4/? 

11 

13 

15 

15 

Example  2 
a  =  v^/30 

/i=l/45,if=l/3.ovlp=5/! 

11 

17 

28 

35 

/i=l/45,if=l/9.ovlp=2/i 

9 

12 

14 

16 

h=\/Ab,H  =  l/lb,o\\p=lh 

7 

8 

10 

10 

/i=l/60,i/=l/10,ovlp=2/2 

9 

11 

14 

15 

h=l/m,H=l/lb.o\\p=lh 

7 

9 

10 

11 

/i=l/60,//=l/20,ovlp=U 

6 

8 

10 

10 

Example  2 
a  =  v/2/100 

/i=l/60.i7'  r  l/15.ovlp=l/j 

8 

11 

18 

20 

/i=l/60.j.  --^  '20.ovlp=U 

7 

9 

13 

16 

This  table  shows  that  the  number  of  iterations  almost  does  not  depend  on 
the  mesh  parameters  as  well  as  the  time  step  size,  especially  for  symmetric 
problems.  For  nonsymmetric,  see  Theorem  3.1.  For  nonsymmetric  prob- 
lems, a  relatively  small  coarse  mesh  size  is  needed  for  reducing  the  total 
number  of  iterations,  that  was  reflected  in  Theorem  3.1  part(3).  The  larger 
the  Reynold's  number  is,  the  smaller  the  coarse  mesh  size  is  needed.  The 
same  as  in  the  stationary  case,  using  less  overlap  can  increase  the  total 
number  of  iterations.  This  can  be  seen  from  the  last  row  of  Example  0. 
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4.2.2      Tests  of  the  modified  additive  Schwarz  method 


e 

1.5 

1.0 

0.5 

0.25 

Example  0 

h=\lAb.H=ll2,.ow\p=bh 

10 

12 

14 

15 

/i=l/60,/f=l/3,ovlp=6/j 

9 

12 

15 

16 

h=l/60,i/=l/10.ovlp=2/j 

12 

24 

34 

35 

;i  =  l/45,i/=l/3,ovlp=/i 

12 

21 

27 

28 

Example  1 

h=\/lbM  =  ll2,.ov\p=2h 

13 

15 

15 

15 

h=l/30,i/=l/3.ovlp=4/i 

13 

14 

16 

16 

/i=l/45.i/=l/3,ovlp=C'. 

13 

15 

16 

17 

/i  =  l/60.if=l/3.ovlp=S- 

12 

15 

17 

17 

/i=l/15.if=l/5.ovIp=: 

14 

18 

19 

20 

/i  =  l/30,i/  =  l/5.ovlp= 

15 

20 

23 

24 

/i=l/45.if=l/5.ovlp= 

13 

19 

24 

25 

h  =  l/60.if  =  l/5.ovlp=-.  . 

13 

18 

23 

25 

e 

1.5 

1.0 

0.75 

0.5 

Example  2 
a  =  \/2/30 

h=l/45.i/=l/3.ovlp=5/i 

13 

22 

31 

42 

h=l/45.i7=l/5.ovlp=3/? 

12 

18 

23 

32 

/i  =  l/45.i/  =  l/9.ovlp=2/i 

10 

14 

19 

30 

/i  =  l/45,if"=l/15.ovlp  =  U 

7 

12 

22 

36 

/i=l/60.i/=l/10,ovlp-2;7 

10 

13 

20 

33 

/j  =  l/60,if=l/15.ovlp=i;i 

8 

13 

24 

41 

Example  2 
a  =  v/2/100 

/i  =  l/45.if=l/9.ovlp=2h 

10 

19 

31 

+50 

/z  =  l/45.if=l/15,ovlp  =  U 

7 

13 

25 

48 

h=l/60,i/=l/15,ovlp=l/? 

7 

13 

25 

+50 

/i  =  l/60.if  =  l/20,ovlp=l/! 

6 

12 

24 

+50 

This  table  shows  that  if  the  factor  r/H'  is  small,  the  results  are  quite 
satisfactory  compared  with  the  corresponding  .4sAf  using  coarse  mesh.  This 
is  ture  especially  for  symmetric  and  mildly  nonsymmetric  problems;  cf. 
Example  0  and  Example  1.  However,  if  t/H'^  is  large,  i.e.  the  time  step  or 
the  number  of  substructures  is  large,  the  number  of  iterations  can  be  large; 
cf.  the  third  row  of  Example  0.  For  problems  with  high  Reynold's  number, 
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the  algorithm  becomes  very  sensitive  to  the  time  step  parameter  e;  cf.  the 
last  column  of  the  second  table.  See  also  Theorem  3.2. 

We  also  note  that  using  less  overlap  can  increase  the  total  number  of 
iterations.  This  can  be  seen  by  comparing  the  first  and  last  rows  in  Example 
0. 

By  comparing  with  t^  -■^S"'^^'  ^^  ^^^  that  this  modified  algorithm  is 
more  sensitive  to  the  m<  arameters  and  the  overlapping  factor. 

4.2.3      Tests  of  the  iterat'-      substructuring  method 


Example  0 

( 

1.5 

1.0 

0.5 

0.25 

;z=l/45.if=l/3 

14 

14 

14 

14 

/i=l/60.if=l/3 

15 

15 

14 

14 

/i=l/60,i/=l/10 

11 

12 

12 

13 

Example  1 

/i=l/45.i/  =  l/5 

IS 

17 

17 

17 

/j  =  l/45.i/  =  l/9 

14 

14 

14 

14 

/r=l/60,i/  =  l/5 

IS 

17 

17 

IS 

/i=l/60,i/=l/10 

13 

14 

14 

14 

Example  2 
(7  =  v/2/30 

/i  =  l/45.ff=l/5 

15 

25 

35 

3S 

h=\l-\.H=\l^ 

13 

IS 

20 

22 

h=l,-ro.H=\/l^ 

10 

12 

13 

14 

h=\/m.H=l/b 

15 

22 

36 

41 

h=l/m.H^l/\Q 

13 

17 

20 

22 

/i=l/60.iy=l/15 

11 

14 

15 

16 

Example  2 
G  =  v^/100 

h=\/AbM=llVo 

11 

15 

20 

22 

/i=l/60,if=l/15 

12 

18 

24 

27 

/i=l/60,i7=l/20 

10 

14 

17 

18 

For  the  symmetric  problems,  we  see  that  the  number  of  iterations  is 
independent  of  the  fine  mesh  size  /i,  but  depends  on  the  coarse  mesh  size  H . 
The  number  of  iterations  is  reduced  when  the  coarse  mesh  size  is  reduced. 
This  is  not  true  for  time  independent  problems.  The  algorithm  is  not  very 
sensitive  to  the  time  step  parameter. 
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For  nonsymmetric  problems,  the  algorithm  becomes  more  sensitive  to 
the  time  step  parameter  and  the  coarse  mesh  size  is  more  important  for 
controUing  the  number  of  iterations. 
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